Background

Over the past few decades there have been numerous publications on the pharmacokinetic (PK) and pharmacodynamic (PD) models of drugs used in anaesthetic practice. These models can predict a drug’s distribution, elimination and clinical effect.

An important use of these models is within target controlled infusion (TCI) systems where the optimal administration of a drug to achieve rapid and stable plasma and effect site concentrations are predicted. The goals of these system are to quickly achieve a clinical effect whilst limiting the side effects of the medication.

Compartment model

The majority of anaesthetic drugs can be described by the compartment pharmacokinetic model with first order elimination. This can take the form of a two or three compartment model with the addition of an effect compartment.

Three

three compartment model

A drug is administered into the central compartment via the intravenous route. The concentration within the central compartment is given by the mass of the drug in the central compartment divided by the volume of the compartment (V1). This is known as the plasma concentration. The plasma concentration rapidly reaches its maximum after the drug administration is stopped. The drug can only be removed from the system via the central compartment though elimination and metabolism.

Once a drug is present in the central compartment it also starts to be distributed to the other two compartments; occuring with rapid and slow equilibration. \(k\) is the equilibrium rate constant and is the rate at which the drug is distributed (units \(min^{-1}\)). The subscript numbers denoting the movement direction; for example \(k_{12}\) is the rate at which the drug is moving from the central compartment to the rapidly equilibriating compartment. \(k_{10}\) is the elimination and metabolism rate from the system.

Once the drug is present in the rapid and slow compartments there is a fraction of drug that moves back into the central compartment. After the administration of a drug into the central compartment stops this movement will maintain a slower decrement of the central compartment concentration. Once in the central compartment this drug can then undergo elimination and metabolism.

Two

The two compartment behaves in the same way as the three compartment model; however, with one less compartment to distribute to.

two compartment model

Some models require the addition of elimination from the second compartment; such as Hofman elimination in cisatracurium modelling. This rate constant is called \(k_{20}\).

Effect site

three compartment model with effect

The central compartment is not where the action of the administered drug occurs. This site of action can be thought of occuring within a seperate effect site which is linked to the central compartment.

The effect compartment is theoretical compartment with neglible volume that links the pharmacokinetic and pharmacodynamics of a drug. The equilibrium between the central and effect site compartment is described by the effect site equilibration rate constant \(k_{e0}\).

\(k_{e0}\) can then be used to calculate the concentration within the effect site compartment. Targetting the effect site concentration rather than the plasma concentration allows for a more rapid attainment of the drug effect at the cost of an overshoot in the plasma concentration.

Mathematical expressions

By applying mathematical models to the compartment model the plasma and effect site concentration after drug administration can be found. Knowing the drug-concentration relationship allows for prediction of dosing required to achieve certain concentrations.

The mathematical models used to create the target controlled infusions are below.

Plasma concentration

The plasma concentration is calculated using the methods outlined by Dubois et al in the PFIM software.

For multiple doses using a three compartment model the plasma concentration is given by:

\[\begin{split}C_p(t) = \left\{ \begin{array}{ll} \sum\limits_{i=0}^{n-1}{\frac{D_i}{T_{inf_i}} \begin{bmatrix} \frac{A}{\alpha}(1-e^{-\alpha T_{inf_i}})e^{-\alpha (t - t_{D_i} - T_{inf_i})} \\ + \frac{B}{\beta}(1-e^{-\beta T_{inf_i}})e^{-\beta (t - t_{D_i} - T_{inf_i})} \\ + \frac{C}{\gamma}(1-e^{-\gamma T_{inf_i}})e^{-\gamma (t - t_{D_i} - T_{inf_i})} \end{bmatrix}} \\ + \frac{D}{T_{inf_n}} \begin{bmatrix} \frac{A}{\alpha}(1 - e^{-\alpha(t - t_{D_n})}) \\ + \frac{B}{\beta}(1 - e^{-\beta(t - t_{D_n})}) \\ + \frac{C}{\gamma}(1 - e^{-\gamma(t - t_{D_n})}) \\ \end{bmatrix} \text{if } t - t_{D_i} \leq T_{inf}, \\ \sum\limits_{i=0}^{n}{\frac{D_i}{T_{inf_i}} \begin{bmatrix} \frac{A}{\alpha}(1-e^{-\alpha T_{inf_i}})e^{-\alpha (t - t_{D_i} - T_{inf_i})} \\ + \frac{B}{\beta}(1-e^{-\beta T_{inf_i}})e^{-\beta (t - t_{D_i} - T_{inf_i})} \\ + \frac{C}{\gamma}(1-e^{-\gamma T_{inf_i}})e^{-\gamma (t - t_{D_i} - T_{inf_i})} \end{bmatrix}} \text{if not.} \end{array} \right\} \\\end{split}\]

where \(C_p(t)\) = plasma concentration at time \(t\), \(D_i\) = total dose of ith infusion, \(T_{inf_i}\) = total duration of ith infusion, \(t\) = current time, \(t_{D_i}\) = start time of ith infusion, \(A\) = phase coefficient of central compartment, \(B\) = phase coefficient of fast compartment, \(C\) = phase coefficient of slow compartment, \(\alpha\) = phase rate constant of central compartment, \(\beta\) = phase rate constant of fast compartment, \(\gamma\) = phase rate constant of slow compartment.

For two compartment modelling the above formula is used with \(C = 0\) so the third (slow) compartment evalutes to 0.

For one compartment modelling both \(B\) and \(C\) are set as 0; with \(\alpha = k_{10}\) and \(A = \frac{1}{v_1}\). \(v_1\) is the volume of the central compartment.

Effect site concentration

The effect site concentration is estimated using the semi-compartment model / direct pharmacodynamic fit detailed by Ki.

Given an estimate of \(k_{e0}\) and observed values of \(C_p\); the following formulas are used recursively starting at time zero when the plasma and effect site concentration are also zero.

Linear solution for an increasing plasma concentration:

\[C_{e_j} = C_{e_{j-1}} \cdot e^{-k_{e0}(t_j-t_{j-1})} + \frac{k_{e0}}{k_{e0} - \lambda_{j_{inc}}} C_{p_{j-1}} \left\{ e^{-\lambda_{j_{inc}}(t_j-t_{j-1})} - e^{-k_{e0}(t_j-t_{j-1})} \right\}\]
\[\lambda_{j_{inc}} = \frac{C_{p_j} - C_{p_{j-1}}}{t_j - t_{j-1}}\]

Log-linear solution for a decreasing plasma concentration:

\[C_{e_j} = C_{e_{j-1}} \cdot e^{-k_{e0}(t_j-t_{j-1})} + (C_{p_{j-1}} - \frac{\lambda_{j_{dec}}}{k_{e0}} ) \left\{ 1 - e^{-k_{e0}(t_j-t_{j-1})} \right\} + \lambda_{j_{dec}}(t_j-t_{j-1})\]
\[\lambda_{j_{dec}} = \frac{\ln{C_{p_j}} - \ln{C_{p_{j-1}}}}{t_j - t_{j-1}}\]

where \(C_{e_j}\) = effect site concentration at time \(t_j\), \(C_{e_{j-1}}\) = effect site concentration at time \(t_{j-1}\), \(t_j\) = time, \(t_{j-1}\) = previous time, \(C_{p_j}\) = plasma concentration at time \(t_j\), \(C_{p_{j-1}}\) = plasma concentration at time \(t_{j-1}\), \(k_{e0}\) = effect compartment equilibrium rate constant

Plasma targetting

In the three compartment model after a bolus injection the plasma concentration is given by:

\[C_{p_{bolus}}(t) = D_{bolus} (A.e^{-\alpha (t - t_D)} + B.e^{-\beta (t - t_D)} + C.e^{-\gamma (t - t_D)})\]

where \(C_{p_{bolus}}(t)\) = plasma concentration at time \(t\), \(D_{bolus}\) = bolus dose, \(t\) = current time, \(t_{D}\) = start time of bolus dose, \(A\) = phase coefficient of central compartment, \(B\) = phase coefficient of fast compartment, \(C\) = phase coefficient of slow compartment, \(\alpha\) = phase rate constant of central compartment, \(\beta\) = phase rate constant of fast compartment, \(\gamma\) = phase rate constant of slow compartment.

An infusion of a drug represents the infinitesimal repetition of smaller bolus doses of \(D\) over time \(t_{inf}\). The plasma concentration of an infusion can be calculated by integrating the above equation.

\[C_{p_{inf}}(t) = D_{{}/{t_{inf}}} \int\limits_{0}^{t_{inf}} {C_{p_{bolus}}(t)}dt\]

where \(C_{p_{inf}}(t)\) = plasma concentration at time \(t\), \(D_{{}/{t_{inf}}}\) = dose per unit time, \(t_{inf}\) = total duration of dose.

Rearranging the equation gives the required dose per unit of time \(D_{{}/{t_{inf}}}\) over time \(t_{inf}\) to reach a plasma concentration \(C_p\) at time \(t\):

\[D_{{}/{t_{inf}}} = \frac{C_p(t)} {\int\limits_{0}^{t_{inf}} {C_{p_{bolus}}}dt}\]

To take into account a situation where other infusions have occured or are currently occuring the required plasma concentration delta is found. This is the difference between the target plasma concentration \(C_{p_{target}}\) and the actual plasma concentration at the specified time to reach the target by \(C_p(t_{inf} + t)\):

\[\Delta C_p = C_{p_{target}} - C_p(t_{inf} + t)\]

where \(C_{p_{target}}\) = target plasma concentration, \(C_p(t_{inf} + t)\) = actual plasma concentration at time \(t_{inf} + t\).

\(C_p(t_{inf} + t)\) is calcuated using the Plasma concentration method.

Therefore the infusion dose per unit time \(D_{{}/{t_{inf}}}\) can be obtained using a target concentration giving \(\Delta C_p\), time of target change \(t\) and duration to achieve the target over \(t_{inf}\):

\[D_{{}/{t_{inf}}} = \frac{\Delta C_p} {\int\limits_{0}^{t_{inf}} {A.e^{-\alpha (t - t_D)} + B.e^{-\beta (t - t_D)} + C.e^{-\gamma (t - t_D)dt}}}\]

Effect site targetting

There are two methods used to target the effect site, original and revised, which are described here.

These are both non-linear problems and are solved by numerical methods using Scipy’s Newton Secant method to find the root of the following functions.

Original

The function takes a plasma concentration target \(C_{p_{target}}\) and effect site target \(C_{e_{target}}\). Then returns the difference between \(C_{e_{target}}\) and the maxima of the effect site concentration curve \(C_{e_{maxima}}\) that results from the initial \(C_{p_{target}}\).

original effect site targetting

By optimising the function by the variable \(C_{p_{target}}\) when the function returns 0 (\(C_{e_{target}} - C_{e_{maxima}} = 0\)) then \(C_{p_{target}}\) is the minimum plasma concentration overshoot required to reach \(C_{e_{target}}\).

Revised

The function takes a plasma concentration target \(C_{p_{target}}\), duration of time to maintain at the target plasma concentration \(t_{inf}\), and effect site target \(C_{e_{target}}\). Then returns the difference between \(C_{e_{target}}\) and the maxima of the effect site concentration curve \(C_{e_{maxima}}\) that results from the initial \(C_{p_{target}}\) over \(t_{inf}\).

revised effect site targetting

By optimising the function by the variable \(t_{inf}\) for a given \(C_{p_{target}}\) when the function returns 0 (\(C_{e_{target}} - C_{e_{maxima}} = 0\)) then \(t_{inf}\) is the duration that the plasma concentration needs to by kept at to reach \(C_{e_{target}}\).

Maintenance infusions

Maintenance infusions are the infusions required to maintain a steady state for a desired plasma or effect site concentration target. In steady state the plasma and effect site should rapidly equilibriate therefore solutions may only need to be performed for the plasma site.

A \(t_{inf}\) is set which is the duration to calculate the maintenance infusion over. The difference between the plasma concentration target \(C_{p_{target}}\) and the actual plasma concentration after \(t_{inf}\) \(C_p(t_{inf} + t)\) is calcuated giving \(\Delta C_p\):

\[\Delta C_p = C_{p_{target}} - C_p(t_{inf} + t)\]

where \(\Delta C_p\) is the plasma concentration difference, \(C_{p_{target}}\) = target plasma concentration, \(C_p(t_{inf} + t)\) = actual plasma concentration at time \(t_{inf} + t\).

Then by rearranging the Plasma concentration equation the required dose per unit time \(D_{{}/{t_{inf}}}\) over duration \(t_{inf}\) to return the plasma concentration back to \(C_{p_{target}}\) can be found:

\[D_{{}/{t_{inf}}} = \frac{\Delta C_P }{ \frac{A}{\alpha}(1 - e^{-\alpha t_{inf}}) + \frac{B}{\beta}(1 - e^{-\beta t_{inf}}) + \frac{C}{\gamma}(1 - e^{-\gamma t_{inf}}) }\]

where \(D_{{}/{t_{inf}}}\) dose per unit time, \(\Delta C_p\) is the plasma concentration to increment by, \(t_{inf}\) is the duration to increment over, \(A\) = phase coefficient of central compartment, \(B\) = phase coefficient of fast compartment, \(C\) = phase coefficient of slow compartment, \(\alpha\) = phase rate constant of central compartment, \(\beta\) = phase rate constant of fast compartment, \(\gamma\) = phase rate constant of slow compartment.

Ke0 ‘tpeak’ method

Traditionally the \(k_{e0}\) was calculated using parametric, nonparametric or sequential PKPD methods. An alternative method is to calculate \(k_{e0}\) using the time to peak effect \(t_{peak}\) as described by Minto et al.

After a patient is given a bolus of a drug \(D_{bol}\) at time zero with the assumption that it immediately reaches the central compartment. The \(t_{peak}\) is obtained from the resulting response-time curve. At \(t_{peak}\) the plasma concentration \(C_p(t_{peak})\) is equal to the effect site concentration \(C_e(t_{peak})\). So that for a three compartment model:

\[C_{e_{bol}}(t_{peak}) = C_{p_{bol}}(t_{peak}) = D_{bol} (A.e^{-\alpha t_{peak}} + B.e^{-\beta t_{peak}} + C.e^{-\gamma t_{peak}})\]

where \(C_{e_{bol}}(t_{peak})\) = effect site concentration at time \(t_{peak}\) after a bolus, \(C_{p_{bol}}(t_{peak})\) = plasma concentration at time \(t_{peak}\) after a bolus, \(D_{bol}\) = bolus dose, \(t_{peak}\) = time to peak effect in seconds, \(A\) = phase coefficient of central compartment, \(B\) = phase coefficient of fast compartment, \(C\) = phase coefficient of slow compartment, \(\alpha\) = phase rate constant of central compartment, \(\beta\) = phase rate constant of fast compartment, \(\gamma\) = phase rate constant of slow compartment.

Substituting in \(k_{e0}\) gives the equation to solve:

\[\begin{split}C_e(t_{peak}) = D_{bol} \begin{bmatrix} \frac{k_{e0}A}{k_{e0} - \alpha} (\alpha e^{-\alpha t_{peak}} - k_{e0} e^{-k_{e0} t_{peak}}) \\ + \frac{k_{e0}B}{k_{e0} - \beta} (\beta e^{-\beta t_{peak}} - k_{e0} e^{-k_{e0} t_{peak}}) \\ + \frac{k_{e0}C}{k_{e0} - \gamma} (\gamma e^{-\gamma t_{peak}} - k_{e0} e^{-k_{e0} t_{peak}}) \end{bmatrix}\end{split}\]

This equation is solved for \(k_{e0}\) using Scipy’s brentq method.