PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Math Comput Model. Author manuscript; available in PMC 2010 October 1.
Published in final edited form as:
Math Comput Model. 2009 October 1; 50(7-8): 959–974.
doi:  10.1016/j.mcm.2009.02.007
PMCID: PMC2768143
NIHMSID: NIHMS147670

A Mathematical Model for the First-Pass Dynamics of Antibiotics Acting on the Cardiovascular System

Abstract

We present a preliminary first-pass dynamic model for delivery of drug compounds to the lungs and heart. We use a compartmental mass balance approach to develop a system of nonlinear differential equations for mass accumulated in the heart as a result of intravenous injection. We discuss sensitivity analysis as well as methodology for minimizing mass in the heart while maximizing mass delivered to the lungs on a first circulatory pass.

Keywords: Compartmental model, circulatory system, drug delivery, differential equations, optimal dosing strategies, sensitivity

1 Introduction

Pharmaceutical companies make significant investment in drug design and development which often primarily entails experiments and clinical trials. Mathematical models for drug uptake and metabolism can, if properly used, be of substantial value in this development process. In this paper we focus first on the process of developing such a model. We then illustrate the use of associated computational methodologies including optimization (maximizing drug delivery to target organs while minimizing amounts at other sensitive sites) and sensitivity analysis (with respect to model parameters and initial conditions. In particular we investigate the potential for first pass effects on the heart that can result from intravenous (i.v.) administration of medicines relative to that for orally administered medicines. The goal of such a model is to determine a strategy for i.v. infusion that will minimize the risk of cardiovascular adverse effects given the chemical properties of the drug. Such models would prove to be very valuable for example in guiding dose administration in the treatment of pneumonia as some anti-bacterial agents carry cardiovascular liabilities. Alternatively, it could guide dose administration for an anti-cancer medicine in the treatment of lung cancer. Orally administered drugs do not pose this risk to patients. Also, model predictions could guide the synthesis of drug candidates. Because the model contains parameters describing drug absorption and metabolism, model predictions could help identify candidate drugs with properties that maximize efficacy in the lungs (i.e., kill bacteria, arrest cancer cells) and minimize the risk of cardiovascular side effects (i.e., tachycardia).

To illustrate the modeling process for drug delivery and metabolism, we consider a class of antibiotics for treatment of pneumonia with the lungs as target organs for delivery. We assume that the drugs are suspected of producing heart damage in organisms being treated. Pre-clinical trials in animals have shown that some antibiotics cause arrhythmias in the hearts of the animals when the drug is administered intravenously, and arrhythmias sometimes occur more frequently with a higher overall dose. This is a serious concern because these arrhythmias pose a heart attack risk. Arrhythmias do not occur when the compound is administered orally, but this may be explained by a lower overall dose.

Although the mechanism behind the drug-related arrhythmias is unknown, one hypothesis is that accumulation of the drug in the cardiac tissue, mostly from the first pass of the drug through the body, is one of the main factors. We develop a preliminary mathematical-pharmacokinetic model to analyze how the drug is distributed throughout the body for different dosage strategies. For a good overview of using differential equations to model components of the cardiovascular system see [5]. The main goal of the model used here is to predict, for a given dosage regimen, the accumulation of the drug in the cardiac and pulmonary tissue from its first pass through the organs. Because we are most concerned with the affects of the drug during its first pass through the system, we will not consider the affects of systemic circulation such as metabolism in the liver. With such models, one should be able to find an optimal dosing strategy that maximizes accumulation of the drug in the lungs, in order to treat pneumonia, while minimizing accumulation in the cardiac tissue.

In the next section we give details for a preliminary first-pass dynamic model for delivery of drug compounds to the lungs and heart. This is followed by simulation results reported in Section 3 along with a description of a computational methodology for choosing optimal injection strategies. A general methodology for sensitivity studies along with sample calculations are presented in Section 4.

2 Model Description

A description of the model compartments, or state variables, along with the equations, parameters, and assumptions of the model will be given here. In compartmental modeling, a compartment can represent a physiological subdivision of a system throughout which the behavior of a substance is uniform [4], and an important part of the development of the model is how the compartments are chosen. A grossly simplified representation of the circulatory system can be seen in Figure 1. After the drug is injected into a vein leading to the vena cavae, it passes through the right atrium, right ventricle, pulmonary artery, lungs, pulmonary vein, left atrium, left ventricle, aorta, coronary artery, and heart tissue.

Figure 1
Simplification of circulatory system.

As can be seen in Table 1, compartments were not created for all of these physiological parts. As the drug passes through the veins and arteries, it is assumed that no amount of the drug is lost, and all of the drug passes through. Inclusion of a vein or artery in the model ensures that the flow rate into the next compartment is represented correctly. Because we are most concerned with the flow rates of the drug into the heart and lung, compartments representing arteries preceding these two organs were included. While one could argue that including compartments for the pulmonary vein and right and left sides of the heart would make the model more accurate, additional compartments would require that more parameters be estimated. In this initial model, we have attempted to balance the trade-off between the number of parameters contained in the model and accuracy of the mathematical representation of the physical system. Another reason the compartments were chosen as depicted is that scientists have the ability to take blood sample measurements at many of the physiological sites included as compartments in our model [19]. These measurements could be used in calibration and validation investigations of the model. A schematic of the model is given in Figure 2, and mathematical descriptions of the rate at which the drug enters/exits the compartments in the model are presented in equation (2.1)equation (2.6). These equations were obtained (under assumptions described below) using mass-balance principles, and their initial conditions are given by (2.7).

Figure 2
Compartment model diagram.
Table 1
Time-dependent state variables.
dMPdt=CU(t)QTVPMP(t)
(2.1)
dMLdt=QTVPMP(t)(aaMT(t)K1+MT(t))ML(t)QTVLML(t)
(2.2)
dMTdt=(aaMT(t)K1+MT(t))ML(t)bMT(t)K2+MT(t)MT(t)
(2.3)
dMAdt=QTVLML(t)QTVAMA(t)
(2.4)
dMCdt=QCVAMA(t)(ccMH(t)K3+MH(t))MC(t)QCVCMC(t)
(2.5)
dMHdt=(ccMH(t)K3+MH(t))MC(t)
(2.6)
MP(t0)=ML(t0)=MT(t0)=MA(t0)=MC(t0)=MH(t0)=0.
(2.7)

The parameter CU (t) in equation (2.1) represents the rate at which the drug is injected at time t, and two different mathematical descriptions of the time-dependent parameter CU (t) for the case of constant injection rate are given in expressions (2.8) and (2.9):

CU(t)={fdif0tτ0ift>τ,
(2.8)

CU(t)={fDTif0tT0ift>T.
(2.9)

In expression (2.8), d represents the dose rate, τ represents the dose delivery time, and f is the fraction of unbound drug in the blood plasma; in expression (2.9), D represents the total dose administered, and T is the dose delivery time. These two alternative descriptions of the injection rate allow for the investigation of different questions regarding dosing regimens. Fixing d and varying τ in (2.8) is analogous to injecting the drug at the same rate for different amounts of time, while fixing D and varying T in (2.9) represents applying the same total dose for different amounts of time.

Sheep are important animal models used for the study of drug kinetics and blood flow, and their body weight and size are similar to that of humans [1]. The values for the parameters in Table 2 represent the circulatory and cardiovascular systems of sheep, and values for many of the parameters were obtained from [19]. (We remark that the specific values obtained from [19] may not be appropriate for use in our first pass model, but can be used to produce simulations illustrating the ideas and methodology presented in this paper.) Of these parameters, only QT and QC were measured directly in efforts reported in [19]; values for VP, VL, and VC were estimated by fitting a model to data under conditions different from those in our investigations here. A value for the volume of the ascending aorta in sheep could not be found; however, a value was found for humans [7]. Using the human ascending aorta volume (.082 L) and a value of 5.0 Liter/min for the total cardiac output for humans [12], we obtained the sheep volume by assuming the ratio of ascending aorta volume to total cardiac output is equal between the two species and then rounded this value to the first significant digit. The parameters that are adjustable and will be interesting in attempts made to examine different dosing strategies include d, τ, D, and T. In Table 2 the parameters are denoted (Bio), (Chem), or (Dose) depending on whether their origin is biological (depending on the organism receiving the drug), chemical (depending on the properties of the drug) or dose related (depending on dosing strategy).

Table 2
Parameter values.

When building a mathematical model, it is important to clearly state all of the assumptions that are made, and a summary of those underlying this model are presented next.

  1. In the model represented by equation (2.1)equation (2.7), each of the compartments have a constant volume and are well-mixed.
  2. After injection, typically into a vein leading to the vena cavae, the drug travels from the injection site to the pulmonary artery, which is the input to the lung, an organ of interest in this work. Enroute to the lung, the distribution of the drug is effected by dilution with other venous flow comprising the total cardiac output QT, dispersion, mixing, and plasma binding. Equation (2.1) accounts for the effective rate of change of the mass of the drug prior to entering the lung. The CU(t) term represents the unmixed rate of injection of the drug. The term VP represents an effective mixing volume that accounts for the dispersion and mixing of the drug enroute to the lung [18]. The term QTVP in equation (2.1) represents the rate of flow out of the pulmonary artery. An important assumption being made in both representations of CU(t) (equation (2.8) and equation (2.9)) is that in the absence of plasma binding (i.e., the fraction of unbound drug in the plasma is 1), the drug enters the pulmonary artery at the same rate it was injected. It is also assumed that the drug injection does not affect blood flow; this assumption is valid if the volume of the solution containing the injected drug is small. With respect to injection of different doses, the current version of the model compares injecting different concentrations of the drug in a small fixed volume of solution (i.e., the total mass injected is allowed to vary but the volume of the solution containing the drug does not change). It is also assumed that neither the vena cavae nor the right side of the heart are affecting the rate at which the drug enters the pulmonary artery.
  3. Similar to the manner in which the effects of the vena cavae and right side of the heart on the drug flow rate were ignored, it is assumed that the pulmonary vein and left side of the heart do not affect the flow of the drug.
  4. The organ compartments corresponding to the heart and lung can each be modeled with two compartments. For each organ, there is a compartment representing the blood supply to the organ as well as a compartment for the organ tissue. When modeling organ compartments in physiology, it is common practice to use separate compartments for the blood and tissue parts [13], and this is the approach we have taken for the two organs in our system.
  5. Equation (2.2) describes the flow into and out of the blood supply of the lung. The term QTVP represents flow into the blood supply of the lung. The exchange between the MP and ML compartments is analogous to an exchange that occurs in a connection between two pipes, and we are assuming here that the flow out of the MP compartment is equal to the flow into the ML compartment. The outflow term in the ML compartment contains the parameter VL representing the volume of the lung, and similar inflow and outflow terms can be found in the MA compartment.
  6. Although the pulmonary artery, lung, and aortic arch receive the entire blood supply of the body, the coronary artery only receives a fraction of this supply which is represented by the parameter QC. The inflow and outflow terms in the MC compartment don’t involve the QT parameter but instead involve the parameter QC which is smaller than QT.
  7. We are left with three nonlinear terms in the model to discuss:
    aaMT(t)K1+MT(t),bMT(t)K2+MT(t),andccMH(t)K3+MH(t),
    which represent absorption into the lung tissue, metabolism in the lung tissue, and absorption into the heart tissue, respectively. The parameters in these terms are heavily dependent upon the particular drug under investigation. The first of these terms is found in equation (2.2) and equation (2.3) and represents the rate at which the drug is absorbed into the lung tissue from the blood supply to the lung. Here a is a parameter representing the maximum absorption rate into the lung with units of (time)−1, and K1 dictates how fast the rate of absorption decreases as the mass of the drug in the lung tissue increases. In the second term, bMT(t)K2+MT(t), b is a parameter that represents the maximum metabolic rate of the drug in the lungs, and K2 determines how fast the metabolic rate saturates. It is known that enzyme-mediated binding and metabolism in the lungs are saturable processes ([10, 18]), and these two terms are represented in Briggs-Haldane Michaelis-Menten form (see [3, 8, 11, 15]). Plots of the absorption and metabolic rates as functions of MT are presented in Figure 3, and values chosen for the parameters are given in the plot. The curve for the absorption rate in the lung intersects the vertical axis at a and then approaches zero as MT goes to infinity. The metabolic rate curve starts out at zero for MT = 0 and then approaches b as MT increases. When the absorption rate curve is greater than the metabolic rate curve, accumulation in the lung will occur; however, when the metabolic rate curve is greater, there will be a net loss of drug mass in the lung tissue. The term (ccMH(t)K3+MH(t)) in equation (2.5) and equation (2.6) represents absorption of the drug into the heart tissue, and here the form was also chosen so that the absorption process is saturable.
    Figure 3
    Absorption and metabolic rates in lung as function of MT for parameter values a = 1, b = 2, K1 = 1, K2 = 1.
  8. While it is known that organs can serve as reservoirs that accumulate drug in their tissue and then release the drug later, this process was omitted in the current version of the model. Because our model is made to study the first-pass effects of the drug through the heart and lung organs, it is reasonable to ignore the release of tissue-bound drug mass.

3 Numerical Results and Optimization Methodology

Typical model outputs MP (t), ML(t), MT (t), MA(t), MC(t), and MH(t) for the compartments are given in Figure 4 for different values of D and T in expression (2.9). In particular, the graphs in Figure 4 depict trajectories for drug concentrations in various model compartments in response to three different sets of dosing parameters D and T. The purpose of this section is to also illustrate the other types of computations that can be performed using our model. These computations can be used to explore questions that may be more difficult or more costly to answer experimentally. There are many goals that can be specified regarding the trajectory of the mass of the drug in the heart and lung. For instance, one may wish to maximize the ratio of peak mass or total mass of the drug that reaches the lung as compared to the heart. Here we will mathematically formulate the goal of, given a fixed dose, finding the best delivery time to maximize these ratios. Note that the total circulation time in humans is about 1 minute [12]. Because our model was developed to analyze the first-pass distribution dynamics of a drug, it is only meaningful to use the model to examine dosage times where the dose will pass from the injection site through the coronary artery in less than one minute (assuming that total circulation time in humans and sheep is similar). It takes time for the drug to move between compartments, and so in order to be conservative, only injection times less than or equal to 30 seconds will be considered here.

Figure 4
Model output for parameter values found in Table 2 with varying values of D and T: D = 100, T = .5 (solid line); D = 50, T = .3 (dashed line); D = 10, T = .1 (dot-dash line).

Expressions (2.8) and (2.9) are only two possible mathematical representations of the drug injection, and they represent injecting the drug at a constant rate. We have freedom to choose the form for the time-dependent parameter CU (t), and other injection strategies are considered here to examine how the input function changes the output of the model. One additional form for the input can be seen in expression (3.1). In a clinical setting, this injection strategy represents constantly increasing the injection rate to a maximum at t = T/2, and then constantly decreasing the rate back down to zero from t = T/2 to t = T. Geometrically, the expression represents an isosceles triangle, and the equation was derived using the formula for the area of a triangle and assuming the maximum height of the triangle occurs at t = T/2. A plot of this input function (along with the constant rate input function and several other considered input functions) can be seen in Figure 5.

CU(t)={4fDTtTif0tT/24fDTtT+4DTifT/2<tT0ifT<t.
(3.1)

Figure 5
Comparison of different injection strategies with T = .5 minutes and D = 100 mg.

Another form for CU(t) is given by expression (3.2); note the addition of two new parameters - namely α, and β. This function was obtained by transforming a beta distribution (which is usually defined between zero and 1) to be defined for t between 0 and T and to have the area under the function be D. Using the transformed beta distribution for the input allows one to examine a large number of shapes for the input function by changing the shape parameters α and β. Visualizations of the transformed beta distribution input for a few different combinations of α and β are also given in Figure 5. Mathematical expressions for two more injection strategies in Figure 5 named the increasing and decreasing injection rate strategies were obtained similar to the manner in which the equation for the triangle injection rate strategy was obtained (equations for these input functions are not written explicitly here). A visualization of how the trajectories for the distribution of the drug in the heart and lung change for the different injection strategies with D fixed at 100 mg and T fixed at 0.5 minutes can be seen in Figure 6.

CU(t)={fDTΓ(α+β)Γ(α)Γ(β)tα1T(1tT)β1if0tT0ifT<t.
(3.2)

Figure 6
Plots of MT and MH for seven different injection strategies shown in Figure 5.

Using well-defined mathematical descriptions of the goal one wishes to achieve with a dosage strategy, one can perform optimizations on cost functions to find the best strategy. Descriptions of some possible cost functionals are given in (3.3)(3.6). When attempting to perform these optimizations, it became clear that the parameter values given in Table 2 cannot be correct. The important dynamics of the distribution of the drug should occur within the total circulation time as long as the injection time is sufficiently short, and this is not the case when using the parameters given in Table 2. To mitigate this difficulty and allow demonstration computations for the optimization ideas with the model, the parameter VL (the volume of the lung) was changed to .25 L and the parameter VC (the volume of the coronary artery) was changed to .01 L; the main purpose of these changes was to cause the drug to move through the system faster.

In (3.3) the objective is to find the dosage time between zero and 15 seconds that minimizes the ratio of peak drug mass in the heart to peak drug mass in the lung; this is equivalent to maximizing the ratio of peak mass in the lung to peak mass in the heart. Using the MATLAB optimization routine fminsearch() with an initial guess of .125 minutes, it was found that shorter injection times minimize the functional corresponding to (3.3). Indeed, as the tolerances of the optimization routine were made more stringent, the optimization routine returned shorter injection times. This suggests that to minimize the ratio of peak masses, the injection interval should be made as short as possible. This is confirmed in Figure 7, where a plot of the ratio of peak masses for different injection strategies is shown. In the plot, one can see that regardless of the injection strategy, the ratio of peak drug mass in the lung to peak drug mass in the heart is maximized when the injection interval is short. However, notice that for longer injection times, some injection strategies perform better than others.

Figure 7
Plot of the ratio of peak drug mass in the lung to peak drug mass in the heart for different dosing strategies and injection times. Injection times were sampled from .01 to .25 minutes with a step-size of .01.

The functional in (3.4) corresponds to finding a dosage time between zero and 15 seconds that minimizes the ratio of total drug load in the heart to total drug load in the lung over one minute (or equivalently, maximizing the inverse). Interestingly, completely different results from those above are obtained when using the functional (3.4). It was found that if the goal is to maximize the total drug load in the lung to total drug load in the heart, longer injection times should be used. When applying the MATLAB optimization routine fmincon() to minimize the ratio of total drug load in the heart to total drug load in the lung with an initial guess of .125 minutes, the optimization routine returns the constraint on the injection time. For example, when the injection time is constrained to be less than .25 minutes as in (3.4), the optimal time returned by the optimization routine is .25 minutes. This suggests that to maximize the ratio of total drug loads in the lung and heart, longer injection intervals do better. Figure 8 shows how the total drug load ratio increases as the injection time increases regardless of the injection strategy. Notice again that the differences between the injection strategies are minimal for small injection times, but the differences grow as the injection time increases. The differences between the results obtained when using the two different cost functionals illustrate the importance of clearly defining the goal of the injection for the distribution of the drug throughout the body.

Min0<T2.5[Maxt(MH(t;T))Maxt(MT(t;T))],
(3.3)

Min0<T.25[01MH(t;T)dt01MT(t;T)dt].
(3.4)

Figure 8
Plot of the ratio of total drug load in the lung to total drug load in the heart over one minute for different dosing strategies and injection times. Injection times were sampled from .01 to .25 minutes with a step-size of .01.

Optimizations performed in (3.5) and (3.6) are two more examples of the types of investigations that can be considered when searching for optimal dosing strategies and further illustrate the importance of clearly stating the goal regarding distribution of the drug during treatment. The cost functional in (3.5) is similar to that in (3.3) in that it corresponds to maximizing the peak mass of the drug in the lung to that in the heart; however, instead of searching for an optimal dosage time, T is constrained to be .25, and the objective is to find optimal shape parameters α and β. An analogous relationship exists between the cost functionals in (3.6) and (3.4) - both represent maximizing the total drug load in the lung as compared to the heart over the first-pass of the drug through the body, but in (3.6) one optimizes the shape of the beta distribution input function for a fixed dosage time. Here the cost functionals do not involve the affects of changing the dose or the dosage time; they only examine changing the shape of the input function. When performing the optimizations in (3.5) and (3.6), α and β are constrained to be greater than or equal to 1, because these values give the input function a unimodal shape for the optimizations

Mina,β1[Maxt(MH(t;T=.25,α,β))Maxt(MT(t;T=.25,α,β))],
(3.5)

Mina,β1[01MH(t;T=.25,α,β)dt01MT(t;T=.25,α,β)dt].
(3.6)

Although different results were obtained for optimal dosage times depending on which cost function was being used, the optimal values for the shape parameters α and β seem to be the same regardless of which ratio is being considered. Optimizations performed on (3.5) and (3.6) revealed that α, β combinations with β equal to one and α greater than one perform best. One can see in Figure 7 that when considering the ratio of peak masses, the beta distribution injection strategy with α = 3 and β = 1 performs better than the others shown. Also, when considering the ratio of total drug loads in the heart and lung, the same values of α and β perform well. A plot of the beta distribution rate with these values of α and β is not shown here, but the input function with these parameters is uniform increasing and convex over the injection interval.

Note that in a clinical setting, it may be difficult or impossible to obtain some of these more complicated shapes for the injection input. However, this methodology was developed so that it may be applied when longer dosage intervals can be considered for a multiple-pass model. It should be emphasized that these results are sensitive to the parameters being used and may change for different parameter values. To obtain realistic predictions, precise parameter values are needed.

The last injection rate strategy examined in this section corresponds to a clinician administering different total amounts of the drug in a series of pulses - i.e., administering the drug at a constant rate, stopping the administration, re-administering the drug at a constant rate and so forth for a given number of times. A plot of this injection strategy can be seen in the top of Figure 9 for 1 pulse and a total dose of 100 mg, 2 pulses (25 mg per pulse) and total dose of 50 mg, and 3 pulses (11.11 mg per pulse) and a total dose of 33.33 mg. These different dosing strategies permit consideration of the effects of administering a dose at a constant rate versus administering half that dose in two equal pulses versus administering 1/3 the original dose in three equal pulses without changing the total time that the dose is being injected. In the lower graph of Figure 9 is a plot of the two different ratios we wish to maximize (peak mass and total load of the drug in the lung as compared to the heart) as a function of the number of pulses. While these simulations suggest that one can force more of a drug to reach the lung relative to the amount that reaches the heart by administering the drug in a series of pulses, better estimates for the parameters in the model are needed before these predictions can be taken seriously. Interestingly, this type of dosing strategy has been investigated by experimentalists, and Seltzer et al., found that they were able to obtain similar therapeutic effects when administering a drug in a bolus or by a series of pulses [16]. Trajectories of the distribution of the drug in the lung and heart over time for different numbers of pulses are plotted in Figure 10. Observe that as the number of pulses increase, the trajectories for both the heart and lung tissue decrease; however, what is important is that the trajectories for the drug in the heart decrease faster causing the two ratios to increase as the number of pulses increases.

Figure 9
Plot of constant rate injection strategy with different number of pulses and value of two ratios (peak mass of drug in lung to peak mass of drug in heart and total drug load in lung to total drug load in heart) as a function of the number of pulses.
Figure 10
MT and MH for constant rate injection strategy with several different pulsing strategies.

4 Sensitivity Analysis

Sensitivity analysis is an important part of model analysis and validation. Sensitivity functions explicitly show the relationship between the parameters and the state variables. If model solutions are very sensitive to a particular parameter, one would want to have a realiable estimate for that parameter, as variations in the value of that parameter will affect the dynamics of solutions. Sensitivity functions are functions of time that describe how a state variable changes with respect to changes in a parameter. To define these sensitivity functions, we will re-state our model in vector notation.

We may define the p-vector of parameters by

θ=(QT,VP,VL,a,K1,b,K2,VA,QC,VC,c,K3,d),
(4.1)

and the n-vector state variable, for a given parameter vector θ and time t, as

x(t;θ)=[MP(t;θ)ML(t;θ)MT(t;θ)MA(t;θ)MC(t;θ)MH(t;θ)],
(4.2)

which allows us to write (2.1)(2.6) in the compact form as

dxdt=g(x(t;θ);θ),
(4.3)

where the vector function g is defined by the right-side of equation (2.1)equation (2.6). This form of the equations will be used in the sensitivity analysis.

The sensitivities of the state variables (MP,ML,MT,MA,MC,MH) to the parameters θ=(QT,VP,VL,a,K1,b,K2,VA,QC,VC,c,K3,d) were calculated. Our sensitivity analysis is described below, using a derivation obtained from [2].

Differentiating with respect to θ on both sides of (4.3), one obtains

θdxdt=θg(x(t;θ);θ),
(4.4)

which implies that the n × p matrix y = [partial differential]x/[partial differential]θ satisfies

ddtxθ=gxxθ+gθ,
(4.5)

or equivalently,

dydt=gxy+gθ.
(4.6)

Numerical solutions of the sensitivity functions y = (yij) = ([partial differential]xi/[partial differential]θj) are calculated by solving (4.3) and (4.6) simultaneously using θ = [theta w/ tilde], where [theta w/ tilde] denotes a given value for the parameter. In our sensitivity analysis we used the parameter estimates listed in Table 2 for [theta w/ tilde]. We note that in order to compute the sensitivity functions the model solutions are required to be continuous and differentiable with respect to the parameters. Because the equation for CU(t) was defined as in (2.8), a step function, we were not able to calculate the sensitivity to the dose time τ. However, with a reformulation of CU(t) as a smooth function of τ, one would be able to obtain sensitivity functions with respect to τ. We have 13 parameters and 6 state variables, hence computing the sensitivity function involves solving 78 first order differential equations (4.6). We can numerically solve this system of sensitivities simultaneously with the solutions of our model.

Relative sensitivities can be found using solutions x(t; θ), [partial differential]x(t; θ)/[partial differential]θ and parameter estimate θ = [theta w/ tilde], by defining

rsij=θ^jxi(t;θ^)xi(t;θ^)θj,fori=1,,n,j=1,,p.
(4.7)

4.1 Sensitivity Results

Figure 11Figure 13 display some of the sensitivity and relative sensitivity solutions we obtained. The sensitivity functions look very different from the relative sensitivity functions. The sensitivity functions show the sensitivity of the state variables to the parameters which includes the effect of the magnitude of the parameter estimates we used. However, the relative sensitivity functions are sensitivities scaled to account for the magnitude of the solutions and parameter estimates. These sensitivity functions change with time, indicating that a state variable may be sensitive to different sets of parameters depending on the time it is examined. This is important because if one is interested in the peak of the state solution, one could investigate to determine the parameters to which the model solutions are most sensitive at the time corresponding to the peak, before the peak, and after the peak. Note that if a sensitivity function for a state variable to a particular parameter was zero for all time, it was not plotted.

Figure 11
Sensitivity and relative sensitivity functions for MT (t; θ) to all parameters.
Figure 13
Sensitivity and relative sensitivity functions for all states to parameter d.

The sign of the sensitivity values corresponds to different effects variations in the parameter will have on the state variables. If the sensitivity function is positive in some time interval, increasing variations in that parameter, in that time interval, will cause the state variable to increase, depending on the magnitude of the sensitivity. Similarly a negative sensitivity value implies that increasing variations in the parameter will cause the state variable to decrease. If the sensitivity function is zero for some time, the state variable is insensitive to that parameter. For example, in Figure 11 (left) we can see that the sensitivity functions for MT (t; θ) are strictly positive over the time span for parameters a, VL, and K1, and strictly negative for VP, b, and QT . Upon close examination, we find that the sensitivity functions for MT (t; θ) are close to zero but strictly positive for K2 and d. Figure 11 and Figure 12 offer a different perspective of our sensitivity analysis results from the results presented in Figure 13. Figure 11 and Figure 12 depict the sensitivity functions for the state variables, MT (t, θ) and MH(t, θ), respectively, to all the parameters. Alternatively, Figure 13 shows the sensitivities of all the state variables to only one parameter, d. These representations contain the same information, just organized differently. The representation of our sensitivity results in Figure 13 is interesting since d is the dose rate, which is a parameter we have freedom to choose as input into our model. We may want to investigate how the variations in d will affect our model solutions over different time intervals. This relates to the analysis of different dosing strategies. For simplicity, we are only considering the constant rate dosing strategy in this sensitivity analysis. From Figure 11Figure 13 we conclude that sensitivity functions contain a lot of information that can be represented in different ways to offer different perspectives in a modeling effort.

Figure 12
Sensitivity and relative sensitivity functions for MH(t; θ) to all parameters.

At any fixed time, sensitivity functions can be ranked by the parameters to which a given state variable are most sensitive. This is essentially what can be seen in Figure 12 and Figure 13. However it is useful to organize this information more concisely in a table. Because the relative sensitivity functions were solved for using a numerical differential equation solver, they are each discretized functions of time:

f=(f(t0),,f(tN))2N+1,

where t0, t1, …, tN correspond to grid points for the numerical solution of the underlying differential equation (4.3) and equation (4.6). In order to rank the relative sensitivities of each state variable to all the parameters, we took the [ell]2-norm of each relative sensitivity function where the [ell]2-norm is defined as |f|2=(i=0N|f(ti)|2)12. The results of these relative sensitivity rankings can be found in Table 3 and Table 4.

Table 3
The [ell]2-norms of relative sensitivities for MP (t; θ), ML(t; θ), MT (t; θ).
Table 4
The [ell]2-norms of relative sensitivities for MA(t; θ), MC(t; θ), MH(t; θ).

From Table 3, we can see that the state variable corresponding to our first compartment, MP (t; θ), is only sensitive (has a non-zero [ell]2-norm) to the parameters that are found in equation (2.1) - the model equation derived from mass balance of that compartment. However, from Table 3 and Table 4, the other state variables are sensitive to the parameters found in their compartmental model equation, and the parameters that appeared in model equations of compartments upstream. The ranking of these parameters does not seem to follow any pattern corresponding to the flow of our compartments.

Sensitivity analysis is very important for compartmental modeling because it clearly shows the dependence of our model solutions on the parameter estimates we used. This sensitivity analysis was for the constant rate dosing strategy defined in (2.8). For a different dosing strategy, new sensitivity functions could readily be calculated.

5 Discussion and Future Work

Observations of the dynamics of the first pass of a drug through the heart led to this model formulation. Compartments were defined, and model equations were obtained using mass balance. Using plausible sheep parameter values found in literature, we made numerical simulations of our model solutions, and analyzed these solutions for various dosing strategies. We also examined the sensitivity of our solutions to the parameter estimates we used. The analysis presented here demonstrates the capabilities of our model, but there is much more work that could be done.

For the purposes of this paper, values for the parameters corresponding to absorption of the drug into the lung (a and K1), absorption into the heart (c and K3), and metabolism in the lung (b and K2) were assigned arbitrary values. To obtain more realistic estimates for these parameters specific to a given antibiotic being developed by a pharmaceutical company, experiments need to be performed. Indicator-dilution studies [14] have been used to quantify uptake of drugs in the lung, and this method could be used to determine more realistic parameter values for absorption into the lung for the drug. More precise values for the other parameters in the model also need to be obtained; many of the parameter values were taken from [19] where the authors obtained these when fitting a model different from ours to data. To obtain more appropriate values for these parameters, the values should be estimated using true first pass experimental data in an appropriate inverse problem formulation. In addition, it is known that the anatomy of the circulatory system has a lot of interpersonal variability [17]. It would be interesting to compare sensitivity analysis predictions with real data on how this variability affects the distribution dynamics of the drug. Most of the parameters in the model are specific to individuals, and clinicians already use knowledge of pharmacokinetic parameters along with the factors that alter them to design optimal dosage regimens [9]. Experimental data on the variability of the pharmacokinetic parameters along with the factors that affect these parameters specific to the drug being developed by a pharmaceutical company would help with making model predictions. Once reasonable model parameters are determined using experimental data for a given drug and organism, one could carry out meaningful simulations to compare drugs that are efficacious and safe; efficacious but with questionable safety; and efficacious but clearly unsafe. In addition, one could carry out simulations aimed at dosing strategies for attaining some minimum concentration in the lung tissue resulting in cure and some maximum concentration in heart tissue resulting in no side effects. That is, one could seek a threshold concentration in lungs such that all higher concentrations are efficacious; and a threshold concentration in the heart such that all lower concentrations are safe. This knowledge may better inform the drug development process. Finally, sensitivity calculations could be performed to determine mechanisms and areas of drug metabolism on which development should focus in future research.

It should be noted that the modeling process is iterative and should involve communication and collaboration with experimentalists. This iterative process is summarized in Figure 14 and discussed in Chapter 1 of [4]. As the next step of the modeling process it would be helpful to obtain experimental data to validate our model (step (vi) from Figure 14). As part of this process, we also need to carefully formulate the statistical model embedded in our mathematical model. Validation of the model using experimental data would allow for comparison of our model with the physiological system we attempted to describe. This most likely will lead to further refinement of our model, and another iteration of the modeling process. Through the communication that takes place between the mathematicians and the experimentalists in the iterative modeling process, the model will more accurately describe the physiological system, and can eventually be used for prediction.

Figure 14
Iterative modeling process.

Acknowledgements

We would like to thank Daniel Beavers (Baylor University), Abhishek Bhattacharya (University of Arizona), Matthew Causley (New Jersey Institute of Technology), and Juanjuan Chai (Indiana University) with whom this effort began as part of a summer research project [6] at the Industrial Mathematical and Statistical Modeling Workshop for Graduate Students sponsored by the Center for Research in Scientific Computation and the Statistical and Applied Mathematical Sciences Institute held at N.C. State University, July 21–29, 2008. Subsequent research was supported in part by the National Institute of Allergy and Infectious Disease under grant 9R01AI071915-05 and in part by the U.S. Air Force Office of Scientific Research under grant AFOSR-FA9550-08-1-0147. KH and NW were supported as CQSB Graduate Fellows during the course of their efforts.

References

1. Adams D, McKinley M. The Sheep. ANZCCART News. 1995;8 Insert.
2. Bai P, Banks HT, Dediu S, Govan AY, Last M, Lloyd AL, Nguyen HK, Olufsen MS, Rempala G, Slenning BD. Stochastic and deterministic models for agricultural production networks. Mathematical Biosciences and Engineering. 2007;4(3):373–402. [PubMed]
3. Banks HT. Lecture Notes in Biomathematics. Vol 6. Heidelberg: Springer-Verlag; 1975. Modeling and Control in the Biomedical Sciences.
4. Banks HT, Tran HT. Mathematical and Experimental Modeling of Physical and Biological Processes. Boca Raton: CRC Press; 2009.
5. Batzel JJ, Kappel F, Schneditz D, Tran HT. Cardiovascular and Respiratory Systems: Modeling, Analysis, and Control. Philadelphia: SIAM; 2007.
6. Beavers D, Bhattacharya A, Causley M, Chai J, Holm K, Wanner N, Wetherington J, Kepler GM, Banks HT, Cintron-Arias A. Fourteenth Industrial Mathematical and Statistical Modeling Workshop for Graduate Students. NCSU: Technical Report CRSC-TR08-21; 2008. Nov, Cardiovascular events associated with oral and IV-administered antibacterial agents.
7. Beneken JEW, DeWit B. A physical approach to hemodynamic aspects of the human cardiovascular system. Physical bases of circulatory transport; regulation and exchange; the proceedings of a conference held in Sept. 1966. 1967:1–45. [PubMed]
8. Briggs GE, Haldane JBS. A note on the kinetics of enzyme action. Biochem. 1925;19:338–339. [PubMed]
9. Brunton LL, Lazo JS, Parker KL. The Pharmacological Basis of Therapeutics. New York: The McGraw-Hill Companies; 2006.
10. Cassidy SS, Scharf SM. Heart-Lung Interactions in Health and Disease. Danbury: Marcel Dekker; 1989.
11. Michaelis L, Menten M. Die Kinetik der Invertinwirkung. Biochem. 1913;49:333–369.
12. Noble A, Johnson R, Thomas A, Bass P. The Cardiovascular System. New York: Churchill Livingstone; 2005.
13. Ottesen JT, Olufsen MS, Larsen JK. Applied Mathematical Models in Human Physiology. Washington D.C: SIAM; 2004.
14. Roerig DL, Kotrly KJ, Dawson CA, Ahlf SB, Gualtieri JF, Kampine JP. First-pass uptake of Verapamil, Diazepam, and Thiopental in the human lung. Anesthesia & Analgesia. 1989;69:461–466. [PubMed]
15. Rubinow SI. Introduction to Mathematical Biology. New York: Wiley-Interscience; 1975.
16. Seltzer JL, Gerson JI, Allen FB. Comparison of the cardiovascular effects of bolus v. incremental administration of thiopentone. British Journal of Anaethesia. 1980;52:527–530. [PubMed]
17. Uflacker R. Atlas of Vascular Anatomy. Second Edition. Philadelphia: Lippincott Williams and Wilkins; 2007.
18. Upton RN, Doolette DJ. Kinetic aspects of drug disposition in the lungs. Clinical and Experimental Pharmacology and Physiology. 1999;26:381–391. [PubMed]
19. Upton RN. A model of the first pass passage of drugs from i.v. injection site to the heart—parameter estimates for lignocaine in the sheep. British Journal of Anaethesia. 1996;77:764–772. [PubMed]
20. Walsh C. Antibiotics: Actions, Origins, Resistance. Washington D.C: ASM Press; 2003.