%% Estimating the Bioavailability of a Drug
% In this example, you will use the parameter estimation capabilities of SimBiology™ to calculate F, the bioavailability, of the drug ondansetron. You will calculate F by fitting a model of absorption and excretion of the drug to experimental data tracking drug concentration over time.
% This example requires Optimization Toolbox™.
%% Background
% Most drugs must be absorbed into the bloodstream in order to become active. An intravenous (IV) administration of a drug is one way to achieve this. However, it is impractical or impossible in many cases.
% When a drug is not given by IV, it follows some other route into the bloodstream, such as absorption through the mucous membranes of the GI tract or mouth. Drugs administered through a route other than IV administration are generally not completely absorbed. Some portion of the drug is directly eliminated and never reaches the bloodstream.
% The percentage of drug absorbed is the bioavailability of the drug. Bioavailability is one of the most important pharmacokinetic properties of a drug. It is useful when calculating safe dosages for non-IV routes of administration. Bioavailability is calculated relative to an IV administration. When administered intravenously, a drug has 100% bioavailability. Other routes of administration tend to reduce the amount of drug that reaches the blood stream.
%% Modeling Bioavailability
% Bioavailability can be modeled using one of several approaches. In this example, you use a model with a GI compartment and a blood plasma compartment. Oral administration is modeled by a dose event in the GI compartment. IV administration is modeled by a dose event in the blood plasma compartment.
% The example models the drug leaving the GI compartment in two ways. The available fraction of the drug is absorbed into the bloodstream. The remainder is directly eliminated. The total rate of elimination, ka, is divided into absorption, ka_Central, and direct elimination, Cl_Oral. The bioavailability, F, connects total elimination with ka_Central and Cl_Oral via two initial assignment rules.
% ka_Central = F*ka
% Cl_Oral = (1-F)*ka
% The drug is eliminated from the Blood_Plasma compartment through first-order kinetics, at a rate determined by the parameter Cl_Central.
%
% Load the project that contains the model m1.

sbioloadproject('Bioavailability.sbproj','m1');

%
%% Format of the Data for Estimating Bioavailability
% You can estimate bioavailability by comparing intrapatient measurements of drug concentration under different dosing conditions. For instance, a patient receives an IV dose on day 1, then receives an oral dose on day 2. On both days, we can measure the blood plasma concentration of the drug over some period of time.
% Such data allow us to estimate the bioavailability, as well as other parameters of the model. Intrapatient time courses were generated for the drug ondansetron, reported in [2] and reproduced in [1].
% Load the data, which is a table.

load('ondansetron_data.mat');

% Convert the data to a groupedData object because the fitting function sbiofit requires it to be a groupedData object.
%

gd = groupedData(ondansetron_data);

% Display the data.
%

gd

% The data have variables for time, drug concentration, grouping information, IV, and oral dose amounts. Group 1 contains the data for the IV time course. Group 2 contains the data for the oral time course. NaN in the Drug column means no measurement was made at that time. NaN in one of the dosing columns means no dose was given through that route at that time.
% Plot the pharmacokinetic profiles of the oral dose and IV administration.
%

plot(gd.Time(gd.Group==1),gd.Drug(gd.Group==1),'Marker','+')
hold on
plot(gd.Time(gd.Group==2),gd.Drug(gd.Group==2),'Marker','x')
legend({'8 mg IV','8 mg Oral'})
xlabel('Time (hour)')
ylabel('Concentration (milligram/liter)')

% Notice there is a lag phase in the oral dose of about an hour while the drug is absorbed from the GI tract into the bloodstream.
%
%% Fitting the Data
% Estimate the following four parameters of the model:
% Total forward rate out of the dose compartment, ka
% Clearance from the Blood_Plasma compartment, clearance
% Volume of the Blood_Plasma compartment
% Bioavailability of the orally administered drug, F
% Set the initial values of these parameters and specify the log transform for all parameters using an estimatedInfo object.

init = [1 1 2 .8];
estimated_parameters = estimatedInfo({'log(ka)','log(clearance)',...
                      'log(Blood_Plasma)','logit(F)'},'InitialValue',init);

% Because ka, clearance, and Blood_Plasma are positive physical quantities, log transforming reflects the underlying physical constraint and generally improves fitting. This example uses a logit transform on F because it is a quantity constrained between 0 and 1. The logit transform takes the interval of 0 to 1 and transforms it by taking the log-odds of F (treating F as a probability). For a few drugs, like theophyline, constraining F between 0 and 1 is inappropriate because oral bioavailability can be greater than 1 for drugs with unusual absorption or metabolism mechanisms.
%
% Next, map the response data to the corresponding model component. In the model, the plasma drug concentration is represented by Blood_Plasma.Drug_Central. The corresponding concentration data is the Drug variable of the groupedData object gd.

responseMap = {'Blood_Plasma.Drug_Central = Drug'};

%
% Create the dose objects required by sbiofit to handle the dosing information. First, create the IV dose targeting Drug_Central and the oral dose targeting Dose_Central.

iv_dose   = sbiodose('IV','TargetName','Drug_Central');
oral_dose = sbiodose('Oral','TargetName','Drug_Oral');

%
% Use these dose objects as template doses to generate an array of dose objects from the dosing data variables IV and Oral.

doses_for_fit = createDoses(gd,{'IV','Oral'},'',[iv_dose, oral_dose]);

%
% Estimate parameters using sbiofit.

opts = optimoptions('lsqnonlin','Display','final');
results = sbiofit(m1, gd,responseMap,estimated_parameters,doses_for_fit,...
                  'lsqnonlin',opts,[],'pooled',true);

%
%% Interpreting Results
% First, check if the fit is successful.

plot(results)

% Overall, the results seem to be a good fit. However, they do not capture a distribution phase over the first hour. It might be possible to improve the fit by adding another compartment, but more data would be required to justify such an increase in model complexity.
%
% When satisfied with the model fit, you can draw conclusions about the estimated parameters. Display the parameters stored in the results object.

results.ParameterEstimates

% The parameter F is the bioavailability. The result indicates that ondansetron has approximately a 64% bioavailability. This estimate in line with the literature reports that oral administration of ondansetron in the 2-24 milligram range has a 60% bioavailability [1,2].
% Blood_Plasma is the volume of distribution. This result is reasonably close to the 160 liter Vd reported for ondansetron [1]. The estimated clearance is 45.4 L/hr.
% ka does not map directly onto a widely reported pharmacokinetic parameter. Consider it from two perspectives. We can say that 64% of the drug is available, and that the available drug has an absorption parameter of 0.4905/hr. Or, we can say that drug clearance from the GI compartment is 0.7402/hr, and 64% of the drug cleared from the GI tract is absorbed into the bloodstream.
%
%% Generalizing This Approach
% lsqnonlin, as well as several other optimization algorithms supported by sbiofit, are local algorithms. Local algorithms are subject to the possibility of finding a result that is not the best result over all possible parameter choices. Because local algorithms do not guarantee convergence to the globally best fit, when fitting PK models, restarting the fit with different initial conditions multiple times is a good practice. Alternatively, sbiofit supports several global methods, such as particle swarm, or genetic algorithm optimization. Verifying that a fit is of sufficient quality is an important step before drawing inferences from the values of the parameters.
% This example uses data that was the mean time course of several patients. When fitting a model with data from more patients, some parameters might be the same between patients, some not. Such requirements introduce the need for hierarchical modeling. You can perform hierarchical modeling can by configuring the CategoryVariableName flag of EstimatedInfo object.
%% References
% Roila, Fausto, and Albano Del Favero. "Ondansetron clinical pharmacokinetics." Clinical Pharmacokinetics 29.2 (1995): 95-109.
% Colthup, P. V., and J. L. Palmer. "The determination in plasma and pharmacokinetics of ondansetron." European Journal of Cancer & Clinical Oncology 25 (1988): S71-4.
% Copyright 2015-2016 The MathWorks, Inc.