www.gusucode.com > simbio 案例源码程序 matlab代码 > simbio/EstimatingTheBioavailabilityOfADrugExample.m

    %% Estimating the Bioavailability of a Drug
% In this example, you will use the parameter estimation capabilities of
% SimBiology(TM), 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(TM).

% Copyright 2015-2016 The MathWorks, Inc.

%% 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
% amoutn 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 adminsitration 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|.

%%
% 
% <<../bioavailability_model.png>>
% 
%%
% Load the project that contains the model |m1|.
sbioloadproject(fullfile(matlabroot,'examples','simbio','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(fullfile(matlabroot,'examples','simbio','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 66% 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 66% 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 66% 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.
%