# Estimate 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

gd = 33x5 groupedData Time Drug Group IV Oral ________ _______ _____ ___ ____ 0 NaN 1 8 NaN 0.024358 69.636 1 NaN NaN 0.087639 58.744 1 NaN NaN 0.15834 49.824 1 NaN NaN 0.38895 44.409 1 NaN NaN 0.78392 40.022 1 NaN NaN 1.3182 34.522 1 NaN NaN 1.8518 28.972 1 NaN NaN 2.4335 25.97 1 NaN NaN 2.9215 22.898 1 NaN NaN 3.41 20.75 1 NaN NaN 3.8744 18.095 1 NaN NaN 4.9668 13.839 1 NaN NaN 5.8962 10.876 1 NaN NaN 7.8717 6.6821 1 NaN NaN 10.01 4.0166 1 NaN NaN 12.08 2.5226 1 NaN NaN 15.284 0.97816 1 NaN NaN 0 NaN 2 NaN 8 0.54951 5.3091 2 NaN NaN 0.82649 14.262 2 NaN NaN 1.0433 19.72 2 NaN NaN 1.4423 21.654 2 NaN NaN 2.0267 22.144 2 NaN NaN 2.5148 19.739 2 NaN NaN 2.9326 17.308 2 NaN NaN 3.3743 15.599 2 NaN NaN 3.9559 13.906 2 NaN NaN 4.9309 10.346 2 NaN NaN 6.1155 7.4489 2 NaN NaN 8.0002 5.1919 2 NaN NaN 10.091 2.9058 2 NaN NaN 12.228 1.6808 2 NaN NaN

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`

compartmentBioavailability 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);

Local minimum possible. lsqnonlin stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.

### 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

`ans=`*4×3 table*
Name Estimate StandardError
________________ ________ _____________
{'ka' } 0.77947 0.1786
{'clearance' } 45.19 2.8674
{'Blood_Plasma'} 138.73 4.5249
{'F' } 0.64455 0.066013

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.