|Home | About | Journals | Submit | Contact Us | Français|
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.
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 . 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.
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 , 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.
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 . 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).
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):
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 . 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 . (We remark that the specific values obtained from  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 ; 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 . Using the human ascending aorta volume (.082 L) and a value of 5.0 Liter/min for the total cardiac output for humans , 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).
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.
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 . 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.
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.
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.
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.
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.
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
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 . 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.
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
and the n-vector state variable, for a given parameter vector θ and time t, as
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 .
Differentiating with respect to θ on both sides of (4.3), one obtains
which implies that the n × p matrix y = x/θ satisfies
Numerical solutions of the sensitivity functions y = (yij) = (xi/θj) are calculated by solving (4.3) and (4.6) simultaneously using θ = , where denotes a given value for the parameter. In our sensitivity analysis we used the parameter estimates listed in Table 2 for . 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; θ), x(t; θ)/θ and parameter estimate θ = , by defining
Figure 11 – Figure 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.
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 11 – Figure 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.
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:
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 2-norm of each relative sensitivity function where the 2-norm is defined as . The results of these relative sensitivity rankings can be found in Table 3 and Table 4.
From Table 3, we can see that the state variable corresponding to our first compartment, MP (t; θ), is only sensitive (has a non-zero 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.
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  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  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 . 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 . 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 . 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.
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  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.