Stability analysis of a diagonally implicit scheme of block backward differentiation formula for stiff pharmacokinetics models

In this paper, we analyze the criteria for the stability of a method suited to the ordinary differential equations models. The relevant proof that the method satisfies the condition of stiff stability is also provided. The aim of this paper is therefore to construct an efficient two-point block method based on backward differentiation formula which is A-stable and converged. The new diagonally implicit scheme is formulated to approximate the solution of the pharmacokinetics models. By implementing the algorithm, the numerical solution to the models is compared with a few existing methods and established stiff solvers. It yields significant advantages when the diagonally implicit method with a lower triangular matrix and identical diagonal elements is considered. The formula is designed in such a way that it permits a maximum of one LU decomposition for each integration stage.


Introduction
Initial value problems (IVPs) governed by the system of ordinary differential equations (ODEs) often arise in the modeling of physical, chemical, and biological systems. It has emerged widely in the areas of medical sciences such as epidemiology, cell physiology, and pharmacology. According to Gear in [1], due to the existence of greatly differing time constants in these systems, a phenomenon known as stiffness is exhibited. Such systems of equations are typically very stable, but often the conventional numerical methods are inefficient due to the severe step length restriction imposed by the numerical stability requirements (see [2]).
In this article, we consider the model of stiff ODEs system in the form of where y T = (y 1 , y 2 , . . . , y s ), β T = (β 1 , β 2 , . . . , β s ) and A is an s × s matrix.
The eigenvalues λ i of the Jacobian matrix J = ∂f ∂ y of the system in (1), evaluated at (x, y), will be used in the most widely heuristic definition of stiffness by Lambert in [3] as follows: The system in (1) is said to be stiff if (i) Re(λ i ) < 0, i = 1, 2, . . . , s, and (ii) max i | Re(λ i )| min i | Re(λ i )| where the ratio [max i | Re(λ i )|] : [min i | Re(λ i )|] is called the stiffness ratio. As stated in [4], a stiff ODEs system is characterized by the property that the ratio of the largest to the smallest eigenvalue is greater than one.
The general solution of (1) takes the form where κ i are the arbitrary constants, c i are the eigenvectors of corresponding eigenvalues λ i , and ψ(x) is a particular integral. Interpreting x to be time, we denote the first term s i=0 κ i e λ i x c i as the transient solution and the remaining ψ(x) as the steady-state solution. Now, let us assume condition (i) is satisfied, which implies that the term s i=0 κ i e λ i x c i → 0 as x → ∞. Let | Re(λ μ )| and | Re(λ υ )| be defined by Re(λ μ ) ≥ Re(λ i ) ≥ Re(λ υ ) , i = 1, 2, . . . , s, so that κ i e λ μ x c i and κ i e λ υ x c i are the fastest and slowest transient, respectively. If the system in (1) is solved numerically and aimed at achieving a steady-state phase, continuing integration is needed until the slowest transient is negligible. Therefore, the smaller | Re(λ υ )| is, the longer the integration period will be. However, for the larger | Re(λ μ )|, a sufficiently small step length is required so that h will lie within the region of absolute stability of the method (refer to [2,3,5]).
According to [6], due to their relative ease of execution, the subclass of diagonally implicit Runge-Kutta (DIRK) methods has become the most commonly used in solving stiff first order ODEs. The fully implicit Runge-Kutta (FIRK) and DIRK methods on Butcher array have the forms as tabulated in Table 1.
The coefficients of a ij in Table 1 constitute a matrix denoted as A. The method with a full coefficient matrix, i.e., FIRK requires to solve a system of (n-dimensional × r-stages) nonlinear equations in each of its integration stages [8]. To reduce the computational cost of evaluating the stages in the FIRK method, [7][8][9][10][11] opted for the DIRK method. As stated in [6], this method is characterized by a lower triangular A-matrix with a ij = 0 for i < j and is sometimes referred to as semi-implicit or semi-explicit Runge-Kutta methods. According to [8], Newton-type iteration is required when the problem is stiff in order to solve the linear systems at each stage with the coefficient matrix written in the form of Table 1 Butcher array for FIRK and DIRK subclasses of the three-stage IRK methods (see [7]) Iha ii ∂f ∂y . It was reported in [8] that if all a ii are equal, the stored LU factorization of such a single matrix can be used repeatedly. This motivates a maximum of one LU decomposition in which the matrix A is written as a product of a lower triangular matrix L and an upper triangular matrix U. It means A is decomposed as A = LU.
Drugs are administered into the body through several routes. The criteria for the selection of a delivery route are the patient's suitability, solubility of the drug, access to the disease's location, and the effectiveness in dealing with the specific disease. According to [12], intravenous, intramuscular, intranasal, intradermal/transdermal, and oral administration are the main drug delivery routes. Among the various delivery routes, oral delivery is the most widely used and commonly employed route of drug administration [12][13][14][15][16]. This is due to the potential advantages of high patient compliance [12,14], ease of administration [12,15], and pain avoidance [16].
To facilitate a meaningful contribution to the pharmacokinetics field, numerous mathematical models for drug delivery systems have been developed, which can be found in [17][18][19][20][21][22][23]. In [21], the established model described the mechanisms of how the human body handles nicardipine (NC)-cyclodextrin complexes injection in the gastrointestinal (GI) tract, distribution in plasma, and metabolism in the liver. Recently, [22] proposed a differential equation model to investigate the dynamic behavior drug resistance incorporating a delay in the process with conditions under which a Hopf bifurcation takes place, which leads to a periodic solution. The models in [23] outline the mechanisms of drug administration in the human body via oral and intravenous routes. To calculate the exact drug concentration in different compartments of blood and tissue medium, the Laplace transformation and eigenvalue method were applied. In this article, we solve the mathematical models of drug diffusion formulated by [20,21,23] numerically. The models were formulated based on the diffusion process by applying Fick's principle and the law of mass action.
The stiffness phenomenon exhibited in the drug diffusion models needed to be treated using a suitable implicit numerical method. Backward differentiation formula (BDF) is the most popular class of implicit linear k-step methods for solving stiff IVPs. To overcome the drawback of the classical BDF method which approximates y n+1 at x n+1 one step at a time, seminal contributions have been made by [24]. The authors employed the block backward differentiation formula (BBDF) to accelerate the integration process. The BBDF produces a block of y n+1 , y n+2 , . . . , y n+k approximations concurrently at each of the algorithms by using earlier blocks with each block containing k points. Some authors have driven the further development of BBDF that proven to improve the performance in efficiency and accuracy upon solving stiff problems (refer to [25][26][27][28]). Although this approach is well established for stiff ODEs, the proof of stability properties has not been extensively studied. Therefore, it is sensible to contribute the relevant proof and modify existing BBDF methods to better accuracy in realizing the newly established scheme's full potential in approximating the solution of stiff pharmacokinetics models.
Dahlquist in [29] introduced the concept of A-stability to address the impact of stiffness on numerically solving IVPs (see [30]). The need for this desirable property, however, imposes a severe constraint on the option of suitable methods for stiff ODEs. This has, therefore, motivated us to design an efficient two-point block method based on BDF in a diagonally implicit manner that holds an A-stable property.
The rest of the paper is outlined as follows. In Sect. 2, we briefly illustrate the derivation of the proposed method. Section 3 provides the stability analysis and some relevant proofs. The mathematical models in which the method is applied are presented in Sect. 4. Section 5 is devoted to the discussion of the numerical results. Finally, Sect. 6 concludes this paper.

ρ-Diagonally implicit block backward differentiation formula
In this section, we derive the two-point ρ-diagonally implicit block backward differentiation formula (ρ-DIBBDF) for the solution of the pharmacokinetics models developed by [20,21,23]. The general form of our method is based on the definition of linear multistep method (LMM) of BBDF as given by Ibrahim et al. in [24] in the form of where m = 2 is order of the method; k = 1, 2 for y n+1 and y n+2 , respectively; ρ is a free parameter; h n+m is step length; α k,k = 1; f n+k = f (x n+k , y n+k ) and β k-1,k = -ρβ k,k .
In the case of the LMM in (2), we associate the linear difference operator L, defined by Expanding y n+j-1 , y n+k and y n+k-1 as Taylor series about x and collecting the terms in (3) gives L y(x n ); h = C 0 y(x) + C 1 h n+m y (x) + · · · + C p h p n+m y (p) (x).
The constants C p are defined as . . .
where s is the step length ratio s = h n+m h n+m-1 . The values of k = 1, 2 in (5) indicate the first and second point, respectively. Without losing the generality, it will be assumed from here on that α 1,1 = 1, α 2,2 = 1, and α 0,2 = 0. Equating the results term by term with (4) yields s(sρ-s-2) 0 , , Next, the coefficients of α j-1,k and β k+m-2,k obtained in (6) are substituted into (2), and by considering h n+m = h n+m-1 for all n, we have s = 1 in (5), which gives the following corrector formula for the two-point ρ-DIBBDF in constant step length:

Stability analysis
In this section, the stability properties of ρ-DIBBDF will be discussed for ρ ∈ (-1, 1). The parameter ρ is restricted to (-1, 1) so that the underlying formula in (2) satisfies the necessary condition for stiff stability. The relevant proof for k = 1 provided in this section is referred to the theorem by [31] (see [32]). The approach adopted to prove some stability properties for k = 2 is similar to that of k = 1.
The following lemma will be used in the next theorem. Proof See [33].
As outlined in [32], a method is said to be A 0 -stable if and only if its corresponding Q(z, μ) is a Hurwitz polynomial for all μ < 0.
Proof Since this theorem satisfies s = 1, i.e., the constant step length formula, then the underlying formula in (7) is A-stable for all ρ in the interval (-1, 1).
To determine whether the numerical method with the chosen ρ is capable of delivering acceptable results, the stability of the method needs to be investigated. Some desirable stability criteria for a numerical method are discussed here, namely zero stability, absolute stability, and convergence. Later, the restrictive requirement on the step length will be studied. For absolute stability purposes, the free parameter ρ is chosen in which it is restricted to (-1, 1) (see [25]). Here, we select ρ = - 3 4 . A detailed discussion for this selection of ρ can be found in [35].

Zero and absolute stability
Denoting hλ byĥ, we have the following definitions pertinent to stiffness. (7) is zero stable if its characteristic polynomial has a simple root at +1 and all the remaining roots lie strictly within the unit circle (refer to [34]). (7) is A-stable if A ⊇ {ĥ| Re(ĥ) < 0}, where A is the region of absolute stability for the method (refer to [5]).
By setting ρ = - 3 4 ,ĥ = 0 and solving (10) with respect to r, it yields the roots, r = 0.01122 and r = 1. Thus, by Definition 4, we conclude that the method is zero stable.
Next, the boundary of the stability region is determined. The boundary locus technique is the most convenient method for finding the region of absolute stability. The region of absolute stability for a range of θ ∈ [0, 2π] is obtained by setting r = e iθ for which |r| ≤ 1. As depicted in Fig. 1, the stable region lies outside the closed contour of the graph. Thus, by Definition 5, the method is considered A-stable.

Convergence
The basic feature that requires an acceptable LMM is that the numerical solution converges to the exact solution y(x) as h approaches zero. The following theorem, the proof of which can be found in [36], sets out the conditions on f (x, y) which guarantee the existence of a unique solution of the IVP in (1) (as in [3]).
Theorem 3 Let f (x, y) be defined and continuous for all points (x, y) in the region D defined by a ≤ x ≤ b, -∞ < y < ∞, a and b finite, and let there exist a constant L as a Lipschitz constant such that, for every x, y, y * such that (x, y) and (x, y * ) are both in D, Then, if y(a) = β is any given number, there exists a unique solution y(x) of the IVP in (1), where y(x) is continuous and differentiable for all (x, y) in D (refer to [3]).
Consequently, this leads to the definition of convergence as given below.

Definition 7
The LMM in (2) is said to be convergent if, for all IVPs in (1) subject to the hypotheses of Theorem 3, we have that lim h→0 y n = y(x n ), holds for all x ∈ [a, b] and for all solutions {y n }, where nh = xa (refer to [3]).
From (7), we have the formula for the two-point ρ-DIBBDF written for the approximate solutions as follows: hf (x n+1 , y n+1 ) + 6 11 hf (x n+2 , y n+2 ), and for the exact solutions as follows: Given that By denoting Y n+ry n+r = d n+r , where r = -1, 0, 1, 2, then (16) where Q = max a≤x≤b |Y (3) (x)|. As h approaches zero, then |d n+k | ≤ |d n | which implies that Y n+k -Y n ≤ y n+ky n for k = 1, 2. Since the condition in (15) is satisfied, we conclude that the proposed method converges.

Step length restriction
To deal with the stability phenomenon as mentioned previously in Sect. 1, we discuss the following definitions given by [5]. The step length restriction is discussed in a similar way as in [37]. (2) is said to be absolutely stable for a givenĥ if all roots of π(r,ĥ) satisfy |r t | < 1, t = 1, 2, . . . , k; otherwise, the method is said to be absolutely unstable.

Definition 9
The LMM in (2) is said to have a region of absolute stability A , where A is a region of the complexĥ-plane, if it is absolutely stable for allĥ ∈ A . The intersection of A with the real axis is called the interval of absolute stability.
From Fig. 1, one could observe that ρ-DIBBDF is absolutely stable except whenĥ ∈ (0, 15.333). This interval is known as the interval of unstable region. As described in Definitions 8 and 9, the A region is defined by the requirement that all the roots of SP(ĥ) have a module of less than 1 for allĥ ∈ A . For the method to be stable, h must lie within a certain range. Thus, we are motivated to find the value of hλ so that |r| < 1. By substituting the endpoint of the interval into SP(ĥ) in (10), we obtain ϕ = 52.85193134ĥ 2 -53.01621516ĥ + 0.1636333333.
On solving (18), we choose the value of |ĥ| < 1. Therefore, the step length is restricted to |H| < 0.003096032803 which is equivalent to Let us consider a stiff problem with λ = -100, then we will have h < 0.003096032803 -100 , which indicates that absolute stability is achieved if the step length chosen is h < 3.096032803 × 10 -5 .

Mathematical model
In this section, we approximate the numerical solutions for five examples of pharmacokinetics models from various literature. A principle based on the law of conservation of mass is applied in the formulation of the mathematical models. It is known as balance law and can be discovered in [38]. Let c(t) denote the concentration of drug in the compartment at time t (hour), hence the balance law applied is dc dt = input rate of drug -output rate of drug, where dc dt is the rate of change for c(t).

Model A
In 1970, in the works of [17], a linear, two-compartment, and open model for drug distribution was presented. Then, a two-compartment pharmacokinetic model was developed by [19] to model the flow of drugs through the body compartments: the GI tract and the circulatory system. The drug is taken orally on a regular basis, resulting in a dosage pulse delivered to one compartment (GI tract). From that point, the drug moves into another compartment (bloodstream) at a rate relative to its concentration in the first compartment. Finally, the drug is metabolized and removed from the blood at a rate corresponding to its concentration there (refer to [20]). Consider the following model consisting of the concentration of drug in the GI tract c 1 (t) and its concentration in the bloodstream c 2 (t) as a function of time t: where D is the drug intake regimen. Noting that a = 2 ln(2) and b = ln(2) 5 , the solution of differential equations in (20) for 0 ≤ t ≤ 6 are given by

Model B
The formulation of drugs with different drug carrier materials is applied widely in the pharmaceutical sciences field to control drug delivery and improve drug release. The authors in [21] had formulated the NC delivery with cyclodextrin (CD) complexes as excipients. The formulation of the NC with hydroxypropyl-β-cyclodextrin (NC/HPβCD) and triacetyl-β-cyclodextrin (NC/TAβCD) modeled in Fig. 2, where k GI and k Plasma are the rate constants in the GI tract and plasma, respectively.
The rate equations corresponding to the scheme in Fig. 2 are where C A and C B are the concentrations in the GI tract and plasma, respectively, and t ∈ [0, 25]. On solving (22), we have The refined rate constants calculated from the predicted (GI) and experimental (Plasma) are listed in Table 2 as follows.

Model C
In this subsection, three mathematical models for drug administration through oral and intravenous routes in the human body are presented. The models have been formulated in the form of IVPs by [23]. The values of the rate constants for the diffusion models are taken from [23] and the eigenvalues λ computed on Maple are listed in Table 3.

Model C(i)
The first model included two compartments of the GI tract and bloodstream, which illustrated the oral drug administration shown in Fig. 3.   The two-compartment mathematical model in Fig. 3 is presented in a linear ODEs approach as follows: The drug delivery model in (1) is presented in terms of its concentration in two compartments: the first compartment in the GI tract c g (t) and the second compartment in blood c b (t). c 0 denotes the initial concentration of drug dosage, while k 1 and k c indicate the inter-compartment rate constant and clearance constant, respectively, where k 1 , k c > 0. By solving equation (24), we have

Model C(ii)
For this intravenous drug model, two compartments are involved viz. blood and tissue as the main exchangers. The drug is injected into the bloodstream and transmitted to the tissue where the drug has a therapeutic effect. A schematic presentation of the intravenous drug administration based on reversible and irreversible rate constants is shown in Fig. 4.
The mathematical formulation of this two-compartment model can be expressed by the following ODEs: in which k b and k t are the rate constants for the blood and tissue, respectively, with c 0 referring to the initial drug dosage and k c denoting the clearance rate. By applying Laplace transformation (see [23]), the resulting differential equations of (26) are solved, yielding where -ξ i are real eigenvalues and are given by

Model C(iii)
The delivery mechanism of Model C(iii) is a one-directional model which is presented in terms of drug concentrations in the following compartments: arterial blood c ab (t); tissue c t (t); and venous blood c vb (t). The drug administration through the three compartments has the schematic view as illustrated in Fig. 5. The scheme modeled in Fig. 5 is a pictorial presentation of the following system of ODEs: The flow rate of the drug from the arterial blood to the tissue compartment is represented by k ab and from the tissue compartment to the venous blood by k t , while the rate of elimination is denoted by k c . As presented in [23], the analytical solution of (28) is obtained as follows:

Numerical simulation
The applicability of our method is then tested on mathematical models established in Sect. 4 to obtain the approximate values. For the interval of t as given in Sect. 4, the numerical results that present the maximum error MAXE and the execution time in seconds TIME for all models are given in Tables 4-10 with  Figure 6 Graph for log 10 MAXE against log 10 TIME of Model A      Tables 4-10 are described as follows: h step length ρ-DIBBDF ρ-diagonally implicit block backward differentiation formula NDIBBDF new diagonally implicit block backward differentiation formula of order 2 (see [27]) BBDF-α block backward differentiation formula-α with the same back values as our method (see [28]) ode15s quasi-constant step size implementation of the numerical differentiation formula ode23s modified implicit Rosenbrock methods of order 2 The numerical results of ρ-DIBBDF compared to NDIBBDF, BBDF-α, Matlab solvers; ode15s and ode23s are displayed in Tables 4-10 for Models A, B, and C. Our results  Figure 10 Graph for log 10 MAXE against log 10 TIME of Model C(i) demonstrate that for step lengths smaller than 10 -2 in Tables 8, 9, and 10, our method gives better performance in terms of efficiency and accuracy compared to other methods. However, for Model C, the numerical results produced slightly larger MAXE when h = 10 -2 . As discussed in Sect. 3, this is because absolute stability will only be achieved if hλ < 0.003096032803. Since the largest λ for Model C(i) is |λ| = 0.9776, then h is required to be less than 3.1669 × 10 -3 . A similar constraint also happened in Models C(ii) and C(iii), where h is restricted to 2.0935 × 10 -3 and 3.1669 × 10 -3 , respectively. It is important to highlight that it is the demand for accuracy and not for linear stability that restricts the step length.
To observe the performance of our method, graphs of log 10 MAXE against log 10 TIME for Models A, B, and C are depicted individually in Figs. 6-12. The performance graphs indicate that the proposed scheme has a uniform pattern of error scaled, which is a significant improvement compared with existing numerical methods. In terms of the execution  Figure 11 Graph for log 10 MAXE against log 10 TIME of Model C(ii) time, Figs. 6-12 display that the ρ-DIBBDF method is capable of performing faster than the comparing methods for all models.

Conclusion
This study has presented a two-point block method based on BDF in a diagonally implicit structure that holds an A-stable property. Stability analysis showed that the proposed method is zero stable, absolutely stable, and convergent. The relevant proof that the method satisfies the condition of stiff stability was also provided. The results provide a basis for the effect of the increase in error due to the step length that operates outside the stability region. In conclusion, the suggested method appeared to perform better than the existing methods of NDIBBDF (formulated in a diagonally implicit manner), BBDF-α (implemented in a fully implicit manner), and the Matlab solvers. Hence, ρ-DIBBDF is significant to serve as an efficient numerical method and as an alternative solver for the 1.93740 ×10 -3 5.50000 ×10 -1 Figure 12 Graph for log 10 MAXE against log 10 TIME of Model C(iii) pharmacokinetics model of ODEs. For future research, ρ-DIBBDF may also be applied to other models arising in the pharmacokinetics field.