Estimate the Bioavailability of a Drug

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:

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