Let us consider the following simple reversible reaction scheme

Where the first step is reversible with a forward rate constant k₁ and reverse rate constant k₂, and the second step is irreversible with the rate constant k₃. The rate constant has the units of the amount per unit time, which are omitted for simplicity in this tutorial. The rate equations are:

Linear ordinary differential equations can be solved using Laplace transformation or the matrix method.[23712] We used the latter to obtain an exact solution because of its simplicity. The above expressions can be written in matrix form:

By setting

$\left|\begin{array}{ccc}-k_1-\mathrm{\lambda }& k_2& 0\\ k_1& -\left(k_2+k_3\right)-\mathrm{\lambda }& 0\\ 0& k_3& -\mathrm{\lambda }\end{array}\right|=0$, where λ is the eigenvalue, and solving for λ, we get the following characteristic polynomial:

Where p = k₁ + k₂ + k₃ and q = k₁k₃. The solutions of the quadratic equation in Eq. (5) are given by

and

The resulting eigenvalues are 0, −α, and −β. These eigenvalues are substituted into the eigenvalue equation to yield the corresponding eigenvectors of

$\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],\hspace{1em}\left[\begin{array}{c}\mathrm{\beta }-k_1\\ k_1\\ \mathrm{\beta }\end{array}\right],\mathrm{\hspace{1em}}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{\hspace{1em}}\left[\begin{array}{c}\mathrm{\alpha }-k_1\\ k_1\\ \mathrm{\alpha }\end{array}\right]$, spectively.

Because the matrix M has three distinct real eigenvalues and their corresponding eigenvectors, we can express the solution as

At t = 0, A = A₀ and B = C = 0, and A + B + C = A₀ at all times. Using these boundary conditions, we can solve Eq. (6) to get C₀ = A₀, C₁ = −C2 and C₁ = A₀/(β – α). Then, the final expressions for the above differential equations (1–3) are

We can compare these results to those for a consecutive two-step reaction, as described in the previous tutorial.[13] We can easily show that the latter has the eigenvalues of 0, −k₁, and −k₃, by setting k₂ = 0 in Eq. (5). The solutions for B and C, equations 8 and 9 in this study, are very similar to those for the consecutive reaction (Equations 7 and 9 in Reference,[13] respectively). Thus, the expressions for the maximum amount of B (B_max) and the time to reach the peak (t_max) can be obtained in a similar way:

and

respectively.

Figure 1 shows representative amount versus time curves of A, B, and C, plotted on a linear and semi-logarithmic scale for three different cases: (a) k₁ = 10, k₂ = 0.1, and k₃ = 1; (b) k₁ = 5, k₂ = 4, and k₃ = 2; (c) k₁ = 1, k₂ = 5, and k₃ = 0.2. When k₁ and k₃ are relatively larger than k₂, that is, the reverse reaction is negligible, linear plots of the amount versus time curves of A and B for the reversible reaction (Fig. 1a) look very similar to those for a consecutive irreversible reaction (see Figure 2 in Reference[13]). The semi-logarithmic plots for the former case (Fig. 1b), however, look different from those for the latter case (see Figure 3a in Reference[13]). In the former case, two curves are always parallel to each other in the elimination phase, whereas they are not in the latter case. The difference can be explained as follows. For an irreversible reaction, the amount of A is described by one exponential term, whereas for the reversible reaction, it is described by two exponential terms. For the latter, the first term in Eq. (7) or (8) decreases faster than the second term over time, because α is always larger than β. This means that the elimination phase is governed by the second term. Therefore, the first term can be neglected when t is sufficiently large. In semi-logarithmic plots of A and B, thus, both slopes for the elimination phase are linear and identical and given by –β/2.303. The y-intercepts for A and B are $\frac{\mathrm{\alpha }-k_1}{\mathrm{\alpha }-\mathrm{\beta }}{A}_0$ and $\frac{k_1}{\mathrm{\alpha }-\mathrm{\beta }}{A}_0$, respectively. When 2k₁ = 2k₂ + k₃, the intercepts are the same, as shown in Figure 2. If 2k₁ < 2k₂ + k₃, the curves A and B will not cross each other. As a third case, we set k₁ = 1, k₂ = 5, and k₃ = 0.2, which satisfies this condition. It is noteworthy that the semi-logarithmic plots for this case (Fig. 3b) strongly resemble those for the consecutive irreversible reaction (see Figure 3b in Reference[13]). In the latter case, the semi-logarithmic plots of plasma concentrations versus time of A and B become parallel straight lines at large t only when the formation rate of B is much slower than its elimination rate; that is, the metabolism is formation rate-limited (FRL).[131415] In addition, the intercept for [B] becomes greater than that for [A] when the distribution volume of B, V_B, is much smaller than that of A, V_A (see Figure 3c in Reference[13]). However, we need to be cautious when we apply this feature to the interpretation of the time profiles of the semi-logarithmic plots of a drug and its metabolite because this feature is commonly observed in reversible metabolic reactions. The pharmacokinetic parameters for each case, including t_max and B_max, are summarized in Table 1.

The area under the curve (AUC) is an important pharmacokinetic parameter, representing total drug exposure over time. The AUC for A and B can be obtained by directly integrating Eqs. (7) and (8) from zero to infinity, respectively, and expressed as

and

Thus, the ratio of AUC_A to AUC_B can be simplified to $\frac{k_2+k_3}{k_1}$. The AUC for C can also be calculated by integrating Eq. (9) from zero to t and given by

Using a time integration method,[1416] we can obtain the relationship between AUC and rate constant or clearance, which has the units of volume (V) per unit time (CL = k × V). By separating the variables in Eqs. (7–9) and integrating both sides with respect to their respective variables, we can write the equations as follows:

They should satisfy the following boundary conditions: [B] = [C] = 0 at t = 0; [A] = [B] = 0 and V_C × [C] = A₀ at t = ∞. By substituting k₁ = CL₁/V_A, k₂ = CL₂/V_B, and k₃ = CL₃/V_C, we get

Eq. (19) can be rearranged as

These relationships are helpful in understanding pharmacokinetic features and metabolic pathways.[141718]

To get a good approximate solution for a reversible reaction scheme, approximation approaches have been developed and are widely used.[121920] Depending upon the reactivity of the intermediate species, we can use either of two approximation approaches. One is a steady-state approximation, in which the intermediate species B is very reactive and hence short lived. This approximation is useful for explaining reactions involving radical species. The other is an equilibrium approximation, in which B is quite stable and long-lived because the equilibrium reaction proceeds faster than the product formation reaction; that is, k₁, k₂ >> k₃. For both approaches, the following condition should be satisfied:

In the steady-state approximation, B is extremely unstable and very short-lived and thus is present in most of the time in extremely small amounts. By assuming that B ≈ 0, we can get a new boundary condition: A + C = A₀. From Eq. (22), Eq. (2) can be rearranged to yield:

For the steady-state approximation, this ratio should be very small because the amount of B is much smaller than that of A. Substituting Eq. (23) into Eq. (1) and re-arranging it for A, we have

By separating the variables and integrating the resulting equation, we get

Using the above expression and the new boundary condition, A + C = A₀, we obtain

Now we can solve Eq. (2) for B using Eq. (26), to get

We plot the exact and approximate solutions for the case where k₁ = 0.1, k₂ = 1, and k₃ = 5 (Fig. 4a). Because $\frac{k_1}{k_2+k_3}=\frac{0.1}{1+5}~0.017\ll 1$, the requirement for the steady-state approximation is satisfied. As the reaction begins, the amount of B quickly increases but is still extremely small (B_max = 0.016 at t_max = 0.72) for most of the time, as expected. As shown in Figure 4a, the approximate solutions (dashed lines) agree well with the exact solutions (solid lines).

For the equilibrium approximation, the reaction needs to be in pseudo-equilibrium. In such a case, we can express

Where K_eq is an equilibrium constant. Because

we also have

Because of the above condition, it is not feasible to directly solve Eqs. (1) and (2). Instead, we need to combine the two equations to get

If we define X = A + B = A₀ – C, from Eq. (28), we can express A and B as

Now Eq. (30) can be solved for X to get

Then, we have

In Figure 4b, we show the plots of the exact and approximate solutions for the case where k₁ = 10, k₂ = 1, and k₃ = 0.1. Unlike in Figure 4a, the amount of B is not negligible (B_max = 0.87 at t_max = 0.44). The reaction rapidly reaches quasi-equilibrium, after which the amount of B slowly decreases. As shown in Figure 4b, the approximate solutions (dashed lines) agree well with the exact solutions (solid lines).

The difference between steady-state and equilibrium approximate approaches lies in the range of their applicability. [20] When we introduce an additional constraint, k₂ >> k₁ into the equilibrium approach and k₂ >> k₃ into the steady state approach, both methods will yield the same solution.[20] This relationship can be further explained by considering an approximate expression for β. Using the following binomial theorem

we can obtain the binomial expansion of β, which yields

If k₁ + k₂ >> k₃ or k₂ + k₃ >> k₁, the second and higher order terms in the expansion can be negligible, and the first term can be further simplified to be the same as the exponent in the equilibrium approximation [see Eqs. (33–35)] or in the steady-state approximation [see Eqs. (25–27)], respectively. For a refined approximation, we can differentiate Eq. (23) for A and B with respect to time t to get

By substituting Eqs. (3), (23), and (37) into Eq. (29) and rearranging the resulting equation, we get

We expect a better approximate solution from the above equation compared with that produced by Eq. (24).

The reversible reaction scheme is also used to describe drug distribution after rapid intravenous injection. Drug distribution within the body is often described by a two-compartment model with central and peripheral compartments.[21] In this model, an injected drug will be either eliminated or distributed within the central compartment, representing intravascular space, and the peripheral compartment, representing non-metabolizing body tissues. The distribution phase can be expressed by two linear differential equations similar to Eqs. 1 and 2.[21] The solution for drug concentration in each compartment can be obtained using a similar approach and described by two exponential terms.

To explain the kinetics of an enzyme-catalyzed reaction, a Michaelis-Menten model has been widely used.[222324] When an enzyme and a substrate form a complex, the binding of the substrate to the active site of the enzyme is non-covalent and reversible. If the conversion of the enzyme-substrate complex to a product is the rate-limiting process, the concentration of the complex can be considered to be constant. Then, the equilibrium approximation will give a hyperbolic relationship between the initial reaction rate, v_0, and substrate concentration [S], known as the Michaelis-Menten equation,

Where V_max is the maximum reaction rate achieved when all of the active sites are occupied by substrates. The Michaelis constant, in the original article,[2223] is given by $K_S=\frac{k_2}{k_1}$. When a steady-state approximation is used, the constant is given by $K_m=\frac{k_2+k_3}{k_1}$. [24]