Reduced models of drug delivery in the presence of fast protein binding

1.   Introduction

Several simple models are currently used to describe the time dependent dynamics of drug concentration in blood. They reflect characteristic features of three different modes of drug delivery: (a) drug administration via injection (e.g., using a syringe); (b) oral drug administration (e.g., via taking a pill); (c) continuous intravenous drug delivery (e.g., via catheter). The equations (and their solutions) for the three cases are shown below ^[1]. For simplicity, throughout this paper we assume that all the parameters and variables are already re-scaled, they are non-dimensional and have some characteristic numerical values associated with them.

(a) Injection model of delivery:

(b) Oral administration:

(c) Continuous intravenous delivery:

In the above formulas, $t$ is time variable, $T$ is the duration of time interval of interest, $x(t)$ is time dependent (free) drug concentration in blood stream; $\alpha$ is the rate constant of drug absorption through the vessel walls and drug elimination (e.g., by liver). In Case(a), $x^*$ is the initial concentration of drug immediately after the injection; in Case (b), the term $\gamma\exp(-\beta t)$ describes the process of drug absorption in the intestines and its transfer into the blood flow after oral administration of medicine (parameters $\gamma > 0$ and $\beta > 0$ are assumed to be known, or they could be estimated from the data); in Case (c), $h$ is the rate of drug inflow into blood stream via a catheter.

In Figure 1 the curves representing behavior of drug concentrations after one injection and after one oral drug administration are shown.

In Figure 2 the results corresponding to multiple administrations are presented together with the solution curve for continuous intravenous drug delivery. The goal is to have the drug concentration for the longest possible time belonging to the therapeutic window, i.e., being below toxicity level and above the level of no efficacy ^[1] (indicated in the figure by dashed lines).

The models mentioned above do not take into account drug-protein binding in blood. The blood proteins include, e.g., albumin, globulin, immunoglobulin, prothrombin, and fibrinogen. The drugs that may bind selectively to some of these proteins include barbiturates, benzodiazepines, penicillin, valproate, phenytoin, warfarin phenytoin, digoxin, and numerous others ^[2]. While experimental and clinical evidence work on drug interactions when two or more of such drugs are administered simultaneously is plentiful ^[3], the modeling aspect of this important problem in not well developed.

There are several characteristic time scales associated with the drug administration and distribution in the presence of drug-protein binding: Forward and reverse binding reactions occur on the time scale of micro- to milliseconds ^[4], drug distribution to tissues and organs happens on the scale of seconds to minutes, and drug action is extended to hours and sometimes days. This means that, since for practical clinical purposes we are usually interested in the processes related to drug distribution and action occurring on the time scale of hours, the fast time scales associated with protein binding may be described by introducing a small parameter and using perturbation methods ^[5] to analyze corresponding model systems of differential equations describing time dependent drug concentrations. Proteins and drug molecules participating in fast forward and reverse binding reactions in blood will be in quasi-equilibrium and the so-called quasi-steady state approximation ^[6,7,8,9] may be applied for derivation of corresponding reduced models. One classical example of this approach is a well-know Michaelis–Menten–Henri approximation ^[10,11,12] widely used in numerous applications including neuroscience and pharmacokinetics. In this paper the boundary function method algorithm is used for the analysis as the most appropriate for the applications which involve exponential behavior within narrow boundary layers ^[13] occurring in the vicinity of the initial instant of time.

In addition to detailed description of drug concentration dynamics via asymptotically reduced models, another important question being addressed here may be formulated as follows. Comparing the graphs in Figure 2, we see that intravenous drug delivery is the most reliable mode of drug delivery which guaranties, under correct dosage defined by the rate of drug infusion into the blood steam via catheter and the knowledge of drug elimination rate constant, the required steady concentration of drug (in the blood stream) belonging to the therapeutic window. Maintaining constant drug concentration in the blood translates into steady supply of drug to tissues and organs to be treated. If the parameters $\alpha$ and $h$ in the model of Case (c) are known, and thus the steady state "saturation" drug concentration $h/\alpha$ is also know, how does parameter $\alpha$ may need to be adjusted to maintain optimal drug concentration level when drug-protein binding has to be taken into account? Also, it is of interest how the rates of delivery of, e.g., two drugs binding competitively to the same protein must be adjusted for the "saturation" concentration of each drug to belong to the corresponding therapeutic window. The main conclusion following from the analysis is that in the presence of fast protein drug-binding, the rates at which the "saturation" concentrations for different drugs are reached will change, sometimes substantially, but the actual "saturation" concentration levels will stay the same for intravenous drug delivery. So, two or more drugs may be administered intravenously using dosages estimated for each drug as if other drugs were absent. This statement refers only to the effect of protein binding on drug delivery and distribution process. When considering possible drug interference, other considerations, like purely chemical drug interactions, may also need to be taken into account.

2.   Materials, methods and results

2.1.   Model reduction in the case where one drug reversibly binds to protein

We start the discussion of methods and results and, in particular, of the boundary function method algorithm ^[13] with a problem of intravenous delivery of just one drug that may bind to a protein. The basic model for Case (c) from the Introduction is supplied with the "generic" forward and reverse binding kinetics scheme ^[14]:

Here we use notations $X$, $Z$ and $Y$ for free drug molecules, unbound (free) protein and bound protein (or, similarly, bound drug), respectively; $k_1$ and $k_2$ are rate constants of forward and reverse binding reactions. Below we will always use small letters, e.g., $x$, $z$, $y$ to represent concentrations of compounds denoted by capital letters, e.g., $X$, $Z$, $Y$, respectively. For continuous intravenous drug delivery the initial concentration of drug in blood stream is $x(0) = 0$; the initial concentration of free protein is $z(0) = z^*$, and initial bound protein concentration is $y(0) = 0$.

In the presence of protein binding, using the Law of Mass Action ^[7], the original model (1.3) is converted into a system of equations for concentrations of species, which may be written as follows:

Corresponding initial conditions are

From the last two equations of (2.2) it follows immediately that

which means that the number of protein molecules, both free and bound, is conserved.

Substituting (2.3) into (2.2), we obtain:

The "slow" characteristic time in the system is associate with the rate of drug delivery $h$ and the metabolism rate constant $\alpha$. The "fast" characteristic time is associated with forward and reverse binding, i.e., numerical values associated with the re-scaled non-dimensional binding - unbinding rate constants satisfy $k_1\gg 1$ and $k_2\gg 1$, which in turn may be expressed using appropriately introduced small parameter $0 < \varepsilon \ll 1$. In particular, we can use the change of parameters

In the above formula and later in the text, notation $\sigma(\varepsilon) = O(\varepsilon^n), n = 0, 1, 2, \ldots$, means that $\lim_{\varepsilon \rightarrow 0}(\sigma(\varepsilon)/\varepsilon^n) = \mbox{const}$.

Using (2.5) the problem (2.4) may be re-written as follows (here we omit tildes to simply notation):

and the requirement $T = O(1)$ means that the non-dimesionalized time interval over which the analysis is going to be performed does not depend on $\varepsilon$.

Problem (2.6) is a standard singularly perturbed system in, so-called, critical case to which the boundary function method algorithm ^[13] may be directly applied. According to the algorithm, the uniform asymptotic approximation to the solution of (2.6) on the closed time interval $0 \leq t\leq T$ may be constructed in the form:

where $\tau = t/\varepsilon$ is a stretched time variable; $\bar{x}_0(t)$, $\bar{x}_1(t) = O(1)$ and $\bar{y}_0(t)$, $\bar{y}_1(t) = O(1)$ are the regular functions in the leading and the first order approximations; $\Pi_0x(\tau)$, $\Pi_0y(\tau)$, $\Pi_1x(\tau)$, $\Pi_1y(\tau)$ are the boundary layer functions of the leading and the first order. Substituting (2.7) into (2.6) and equating separately regular and boundary layer functions multiplying like powers of $\varepsilon$, we arrive at the problems for the terms in (2.7). Characteristic features of boundary layer functions: (ⅰ) They are needed to compensate for the discrepancies introduced by the regular functions in the initial conditions; (ⅱ) These functions must decay to zero as the stretched variable $\tau$ tends to infinity (this condition allows one to extend the asymptotic method, working naturally for the linear differential equations, to the case of nonlinear differential equations and systems). For problem (2.6), as well as for a more complex problem considered in the next section, all the conditions needed to justify the asymptotics and estimate the remainder terms in (2.7), i.e., the asymptotic order of an error defined as a maximum of an absolute value of the difference between the constructed asymptotics and the "exact solution" taken over the time interval of interest, are satisfied (see ^[13]). In what follows, we will limit our analysis to the asymptotic terms of the order $O(1)$ only.

In the leading order approximation, we get the following set of problems for the terms of asymptotic expansion. For regular functions $\bar{x}_0(t)$ and $\bar{y}_0(t)$:

and thus,

Equation (2.8) is obtained just by setting $\varepsilon = 0$ in (2.6), and (2.9) is the quasi-equilibrium relationship (quasi-steady state between the drug and the protein molecules which exists due to the presence of fast forward and reverse binding), showing that $\bar{x}_0(t)$ and $\bar{y}_0(t)$ are dependent on each other. To get the differential equation for free drug concentration approximation in the leading order, $\bar{x}_0(t)$, we need to discuss the problem for regular functions in the first order approximation:

The problem (2.10) is just a non-homogeneous system of linear algebraic equations for $\bar{x}_1(t)$ and $\bar{y}_1(t)$ whose solvability condition (obtained by just adding the two equations in (2.10)) produces the differential equation involving derivatives of $\bar{x}_0(t)$ and $\bar{y}_0(t)$:

Substituting the expression (2.9) for $\bar{y}_0$ in terms of $\bar{x}_0$ into (2.11) and performing differentiation, we arrive at the differential equation for $\bar{x}_0$:

or

Equation (2.12) must now be solved with the initial condition determined simultaneously with the leading order boundary functions $\Pi_0x(\tau)$ and $\Pi_0y(\tau)$, which satisfy the nonlinear system of equations:

conditions at infinity:

and, together with the leading order regular functions satisfy the initial conditions:

The last equation of (2.14) was transformed using (2.9) and relationship $\Pi_0y(\tau) = -\Pi_0x(\tau)$ following directly from (2.13) and conditions at $\tau = \infty$. The system of two nonlinear algebraic equations of (2.14) for two unknowns has a unique solution: $\bar{x}_0(0) = 0$ and $\Pi_0x(0) = 0$, and thus, $\bar{y}_0(0) = 0$ and $\Pi_0y(0) = 0$. Thus, the problem formulation for $\bar{x}_0(t)$ is now complete: It consists of the differential equation (2.12) and zero initial condition:

The resulting problem for $\Pi_0x(\tau)$,

has only a trivial solution $\Pi_0x\equiv 0$.

Thus, the leading order approximation of the solution of the original problem (2.6) is given by

where $\bar{x}_0(t)$ is the solution of (2.12) with zero initial condition (2.15).

Qualitative behavior of the solution of (2.12) is similar to that of (1.3): For both equations there exists only one stable steady state: $\bar{x}_0 = h/\alpha > 0$ for (2.12) (and $x = h/\alpha > 0$ for (1.3)). The transition to that steady state for the solution of (2.12) is slower compared to that for (1.3) because for the factor multiplying the term $(\alpha \bar{x}_0 -h)$ on the right-hand side of the reduced model equation we have:

Let us emphasize that both the expression for quasi-steady state (2.9) and differential equation (2.12) contain only the ratio of rate constants $(k_2/k_1) = (\tilde{k}_2/\tilde{k}_1)$, which means that if the experimental data are collected on the time interval $[0, T]$, where $T = O(1)$, then the values of these constants cannot be estimated separately from each other, but only their combination $K = k_2/k_1$ can be reliably estimated.

Comparison of the two solutions mentioned above for a particular choice of parameter values is shown in Figure 3.

The presented analysis may be repeated for the injection (Case (a)) and oral administration (Case (b)) modes of drug delivery. The reduced model Eq (2.12) for these cases will stay practically the same with the following changes: The term $(\alpha \bar{x}_0-h)$ on the right-hand side of (2.12) must be substituted for by $\alpha \bar{x}_0$ (in Case (a)) and by $(\alpha \bar{x}_0-\beta\exp(-\gamma t))$ (in Case (b)). The resulting model equation for Case (b) must be solved with initial condition from (1.2). However, the initial condition for model equation in Case (a) will be different from that in (1.1); to find $\bar{x}_0(0)$ for this case we must solve the system similar to (2.14), now containing $x^*$ from (1.1):

The only physiologically meaningful (positive) solution of (2.18) is

This initial condition belongs to the "slow manifold" described by (2.9); the resulting value of $\bar{x}_0(0)$ for free drug is less than $x^*$ since some portion of the originally administered drug is bound to protein at a fast time scale. The fast transition from the original initial condition $x^*$ to $\bar{x}_0(0) < x^*$ within a narrow initial boundary layer is described by the exponentially decaying boundary function $\Pi_0x(\tau)$ obtained as a solution of equation similar to (2.16) with the initial condition $\Pi_0x(0) = x^*-\bar{x}_0(0) > 0$, where $\bar{x}_0(0) > 0$ is given by (2.19).

The comparison of sample concentration curves for a specific choice of parameters is shown in Figure 4. We compare the numerical solutions of the original models without drug-protein binding, (1.1) for Case (a) and (1.2) for Case (b) (solid curves), with the solutions of the corresponding reduced models of type (2.12) obtained in the case of fast forward and reverse binding with protein (dashed curves). Numerical solutions are obtained for a sample set of parameter values. For Case (a): $\alpha = 1$, $z^* = 1$, $k_1 = 100$, $k_2 = 100$, so that $\varepsilon = 0.01$; for Case (b): $\alpha = 2$, $\beta = 1$, $\gamma = 4$, $z^* = 1$, $k_1 = 100$, $k_2 = 100$, so that $\varepsilon = 0.01$.

We note that the effect of protein binding on time dependent drug concentrations in a blood flow is different for different modes of delivery. In the case of injection, as well as in the case of oral administration, the peak value of drug concentration is lower when protein binding is present, which may lead to the optimal dose miscalculation (the drug dose in the presence of drug binding must be increased). Continuous intravenous delivery via catheter is the optimal mode of drug delivery which allows one to keep the required constant drug concentration in blood over long time intervals.

Next, let us apply the methodology outlined in the current subsection to a more complex problem involved with modeling two drugs administered simultaneously via catheter directly into bloodstream (intravenous drug delivery) in the case of competitive binding to the same protein. The reduced model produced in the next subsection is novel: It is not found in the previously published literature. We will present the major features of the analysis omitting some of the cumbersome detail.

2.2.   Model reduction in the case where two drugs bind competitively to the same protein

Consider an extension of the previous model where two drugs can reversibly bind to a protein. For each of the drugs basic model of Case (c) from the Introduction holds, with different values of parameters $\alpha$ and $h$: We will refer to them as $\alpha_1$ and $h_1$ (for the first drug) and $\alpha_2$ and $h_2$ (for the second drug). Each basic model is also supplied with the "generic" forward and reverse binding kinetics reaction scheme involving the same protein molecules:

Here $X_1$ and $X_2$ are the two drugs being administered, $Z$ represents the unbound (free) protein, $Y_1$ and $Y_2$ are the protein molecules bound to the first and the second drug, respectively. The rate constants $k_1$, $k_2$, $k_3$ and $k_4$ of the forward and reverse binding reactions are, in general, different.

Same as in the case of a simpler model described earlier, the system of equations for concentrations of species $x_1(t)$, $x_2(t)$, $z(t)$, $y_1(t)$ and $y_2(t)$ may be written as follows:

The initial conditions reflect the situation where the initial concentrations for both drugs are zero, the initial protein concentration is $z(0) = z^* > 0$ and all protein molecules are unbound:

It follows from the last three equations of (2.21) and initial conditions (2.22) that

which means that the number of protein molecules, both free and bound, is conserved.

Substituting (2.23) into (2.21), we obtain:

In the presence of fast forward and reverse binding, i.e., when the numerical values associated with the re-scaled non-dimensional binding–unbinding rate constants satisfy $k_1\gg 1$, $k_2\gg 1$, $k_3\gg 1$, $k_4\gg 1$, same as before, a small parameter $0 < \varepsilon \ll 1$ may introduced. Then, after applying the change of parameters

the system (2.24) may be re-written in the form of a singularly perturbed problem in the critical case ^[13] (here, same as before, we omit tildes to simply notation):

Asymptotic approximation for the problem (2.26), (2.22) can be constructed using the boundary function method algorithm ^[13] in the form similar to (2.7), but now the regular functions and boundary layer functions must be constructed for four variables: $x_1$, $x_2$, $y_1$ and $y_2$. Here we are only interested in the leading order approximation. Let us use the following notations for the leading order terms:

where $\tau = t/\varepsilon$ is a stretched time variable. Substituting (2.27) into (2.26) and equating separately regular and boundary layer functions multiplying like powers of $\varepsilon$, we arrive at the problems for the terms in (2.27). In the leading order approximation, by setting $\varepsilon = 0$ in (2.26), for $\bar{x}_{10}$, $\bar{x}_{20}$, $\bar{y}_{10}$, $\bar{y}_{20}$, we obtain:

This is a system of two nonlinear algebraic equations from which $\bar{y}_{10}$ and $\bar{y}_{20}$ may be expressed in terms of $\bar{x}_{10}$ and $\bar{x}_{20}$:

The quasi-equilibrium relationships (2.29), which hold due to the fast forward and reverse binding, indicate that in the leading order approximation the drug and bound protein concentrations $\bar{x}_{10}(t)$, $\bar{y}_{10}(t)$ and $\bar{x}_{20}(t)$, $\bar{y}_{20}(t)$ are dependent on each other. The system of differential equations for $\bar{x}_{10}$, $\bar{x}_{20}$ is obtained as follows. Let us add the first and the third equations in (2.26); also let us add the second and the fourth equations in (2.26). These are the equivalent transformations, and so, the resulting two equations must hold for the original unknown functions as well as for their approximations, i.e., for the regular functions in the leading order approximation:

These equations are an analog of the solvability condition (2.11) obtained earlier for the previous simpler example. Substituting the derivatives of $\bar{y}_{10}$ and $\bar{y}_{20}$, computed by differentiating (2.29), into (2.30), we obtain

Let us re-write (2.31) in the form

where for the elements of matrix $A$ we have:

Simple calculation indicates that $\det A > 0$, and so, inverse matrix $A^{-1}$ exists and may be used to produce from (2.32) a system of nonlinear differential equations for $\bar{x}_{10}(t)$, $\bar{x}_{20}(t)$ resolved with respect to derivatives:

The elements of

in terms of $\bar{x}_{10}(t)$, $\bar{x}_{20}(t)$ may be readily written out, but we will not do it here to shorten the presentation.

System (2.33) must be solved with the initial conditions obtained together with the construction of the leading order boundary functions. It can be shown, similar to the case of a simpler example discussed earlier, that all the $\Pi$-functions in the leading order approximation are zero, and that the initial concentrations for the regular terms in the leading order approximation are zero as well:

Qualitative behavior of the solutions of coupled system (2.33) is similar to that of two decoupled equations of type (1.3) where different values of parameters $\alpha$ and $h$ are used for two different drugs. There exists only one stable steady state: $\bar{x}_{10} = h_1/\alpha_1 > 0$, $\bar{x}_{20} = h_2/\alpha_2 > 0$ to which the solution starting at the initial state (2.34) tends as time increases. For each of the two drug concentrations the transition curves are illustrated in Figure 5: (A) We compare the numerical solution of the original model without drug-protein binding consisting of two equations of (1.3) containing different parameters (solid curves) with the solution of the corresponding reduced model (2.33), (2.34) obtained for the case of fast forward and reverse competitive binding of two different drugs to the same protein (dashed curves). All parameters, concentration and time variables are non-dimensional and are assigned some particular numerical values: $\alpha_1 = 2$, $\alpha_2 = 1$, $h_1 = 1$, $h_2 = 2$, $z^* = 1$, $k_1 = 100$, $k_2 = 100$, $k_3 = 100$, $k_4 = 100$, so that $\varepsilon = 0.01$. Random parameter value choices do not have any special physiological significance. (B) Comparison of the numerical solution of the original singularly perturbed model, consisting of system (2.26) supplied with zero initial conditions, and obtained for the numerical values of parameters shown above, with the reduced model solution (leading order approximation for the regular functions obtained by solving (2.33) and (2.34)). The two solutions are practically indistinguishable for the value of $\varepsilon = 0.01$. The graphs exhibit the type of behavior similar to that shown for one drug in Figure 3.

The presented analysis shows that in clinical applications where multiple drugs binding to the same proteins are administered intravenously via catheter, the steady state (saturation) concentrations of different drugs in the blood flow are not affected by presence of other drugs. However, the time characteristics of transition from initial zero drug concentrations to saturation levels are affected by the rate constants of forward and reverse binding reactions, and in the case of fast binding–unbinding the time dependent concentrations are described by (2.12) with zero initial condition (2.15) for the case of one drug, nd by (2.33) with zero initial conditions (2.34) for the case of two drugs.

3.   Discussion and conclusions

The procedure for deriving reduced models of drug delivery in the presence of fast reversible drug–protein binding is presented and illustrated using two examples. One example involves continuous intravenous drug delivery via catheter of one drug binding to a protein. Another, more complex example deals with administering of two drugs competitively binding to the same protein. The derived reduced models allow one to trace time dependent free drug concentrations during the transition periods between the start of administration (with zero initial free drug concentration) and reaching the steady state (with saturation drug concentration in blood flow). The same approach also works for studying other modes of drug delivery, such as administration via injection and oral administration.

The results presented in this paper have both practical and theoretical significance. On the one hand, they may be used for constructing a detailed time dynamics of drug delivery and distribution for particular choices of commonly used drugs. The presented algorithm elucidates an effective general approach to analyzing drug delivery problems for more complex drug-protein forward and reverse binding reactions as well. Derived formulas may be used directly in the studies related to design of new drugs, and modeling their delivery and distribution. On the other hand, the reduced model formulas may be handy for practical clinical applications, guiding the calculation of optimal drug doses for various delivery modes. Correct drug dosing in the presence of drug-protein binding is a long standing problem in pharmacokinetics ^[15]. It is especially pressing in the common clinical situations such as hypoalbuminaemia ^[16], where for a variety of reasons the protein, e.g., albumin, concentration in blood decreases. This is expected to lead to changes in free drug concentration in blood and, in turn, has strong influence on the therapeutic effects and possible side effects (e.g., toxicity for high concentrations) of administered drug. The standard way of analysis for such problems, i.e., determining the correct dosing for various changing protein concentrations in blood, involves the so-called "in vitro" approach, where the drug and protein concentrations are assumed to be static, leading to a steady state fractions of free and bound drug ^[15] after delivery. When the steady state concentration of protein changes, this leads to a new steady state distribution between free and bound fractions of the drug. The tables produced for dose estimates under such artificial conditions are not taking into account different possible modes of drug delivery and they also ignore the dynamic changes of drug concentrations due to different administration modes and drug metabolism. In fact, in the current widely used approach to optimal dosing determination, the situation where the free drug concentration is reduced due to the presence of drug-protein binding when the drug is administered either orally or via single injection, as illustrated in Figure 4, is automatically extended to the intravenous mode of delivery case ^[15], which must not be done. For intravenous drug delivery mode the saturation drug concentration in blood is independent of the observed protein concentration even in the clinical cases characteristic of hypoalbuminaemia (see illustration of this effect in Figures 3 and 5(A)). When the results of the analysis presented in this paper are known, they seem to be very simple and even obvious. However, without this analysis the conclusions are not self-evident, as the current clinical practice shows.

Other asymptotic algorithms, e.g., Matching technique ^[5], may be used for derivation of the presented reduced models as well. The boundary function method, however, is preferable for its simplicity.

Use of Generative-AI tools declaration

The author declares he has not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

No external funding was used. The author would like to thank the anonymous reviewers for their valuable comments and suggestions which led to substantial improvement of the final version of the manuscript.

Conflict of interest

The author declares no conflict of interest.