Rich dynamics of a three-tiered anaerobic food-web in a chemostat with multiple substrate inflow

1.   Introduction

Coronary artery disease (CAD) is a condition where plaque builds up inside the coronary arteries, which are the blood vessels that supply oxygen-rich blood to the heart muscle. As the plaque accumulates, it can narrow or block the arteries, reducing blood flow to the heart and causing chest pain or discomfort, shortness of breath, fatigue, and other symptoms. CAD can also lead to more serious conditions, such as heart attack or heart failure. Treatments for CAD, including angioplasty, vary depending on the severity and extent of the disease. In some cases, a stent may be placed during angioplasty. There are two main types of stents: bare-metal stents and drug-eluting stents (DESs). Bare-metal stents are made of metal and are effective at keeping the artery open, but they can sometimes cause re-narrowing of the artery, called restenosis. DESs are coated with medication that helps prevent re-narrowing and improve long-term outcomes. In order to model the drug delivery from the DES to and through the arterial walls, it is necessary to study the biological structures of the arteries. There are three layers that comprise the arterial wall (see Figure 1), starting from the inside of the wall: intima, media, and adventitia. A thin layer of endothelial cells, called the endothelium, lines the inside of the intima. They are in contact with the blood and control the relaxation and contraction of the artery, as well as prevent the smooth muscle cells in the media from proliferating. The media is composed of smooth muscle cells, collagen, and elastic fibers that help to regulate blood pressure and flow. The smooth muscle cells are the targets for the drug delivery. The adventitia is composed of connective tissue that supports the artery. It is filled with tiny blood vessels called vasa vasorum, which supplies blood to the adventitia and acts as a clearance mechanism for drugs released into the artery wall.

Many multi-layer models have been proposed to study the pharmacokinetics in the arterial wall. Among the one-dimensional models, we first mention the model proposed in ^[2], which consists of a diffusion equation in the drug-coating region and a diffusion-advection-reaction equation in the arterial wall region. The coupling is achieved by applying interface conditions. In ^[3], the authors further took the intracellular concentration into account, along with extending an early model to include the adventitia layer as well. In ^[4], the authors studied the 2-layer model from ^[3] and provided an analytic solution in some special case. This two-layer model is the focus of this paper.

High-dimensional models have been studied as well. We refer the reader to ^{[5,6,7,8,9,10,11,12,13]} for more details. Here, we focus our attention on the model proposed in ^[4] since it models the intracellular concentration separately. The reader is also referred to ^[14] for a review of different models.

The remainder of this paper is organized as follows. In Section 2, we introduce the one-dimensional model and state its weak formulation. In Section 3, we present a complete PDE analysis for the model, which includes derivation of a priori energy estimates and the establishment of its well-posedness by using the Galerkin method with a compactness argument. In Section 4, we present a complete finite-element numerical analysis for the PDE (partial differential equation) model, followed by the numerical results given in Section 5. In Section 6, we first introduce a generalized two-dimensional model and then sketch some PDE and numerical analysis results for the proposed model. Finally, the paper is completed with a few concluding remarks given in Section 7.

2.   Mathematical model and its weak formulation

In this section, we first introduce the one-dimensional drug-release model from ^[4] and then present its weak formulation. We note that several geometric simplifications were adopted when establishing this 1-d model. First, the endothelium is usually severely damaged after the stent insertion; it is therefore omitted. Second, the intima, when devoid of the endothelium, has a structure that is similar to the media and will thus be absorbed into the media region in the model. Third, the adventitia is omitted in the model since research shows that it does not have a large effect on the drug concentration in the media region. See Figure 2 for a schematic diagram.

Let $c$, $c_1$, and $c_2$ denote, respectively, the concentrations of the drug in the stent coating, in the extracellular matrix, and in the smooth muscle cells. The stent concentration is governed by the diffusion equation, the extracellular concentration by the diffusion-advection-reaction equation, and the intracellular concentration by a linear Ordinary Differential Equation (ODE). A no-flux boundary condition is imposed at the lumenal boundary ($x = -l$). The stent and wall concentrations are coupled through the continuity of mass flux, as well as the Kedem-Katchalsky equation at the interface ($x = 0$). The system is then non-dimensionalized. We refer the reader to ^[4] for a detailed explanation. However, we note that the original model proposed in ^[4] imposes a boundedness condition on the solution, whose main purpose is to help one to obtain an analytic solution, but this restriction may not be appropriate from the PDE point of view and, more importantly, it is difficult to approximate numerically. Hence, we chose to replace this boundedness condition by imposing a no-flux boundary condition on the adventitial boundary ($x = 1$), under the assumption that no drug escapes through the adventitial boundary. We note that this is an idealized situation. It is also common to impose the homogeneous Dirichlet boundary condition there, under the assumption that the drug concentration would be negligible at the far end of the arterial wall. Between the two idealized situations, we chose the former, with the understanding that the analysis of the system would not be affected aside from having slightly different solution spaces.

Specifically, let ${\Omega^s}: = (-l, 0)$ and ${\Omega^m}: = (0, 1)$; our one-dimensional drug-release model is given as follows:

where $\partial_t$ denotes the partial derivative in time $t$ and the sub-index $x$ represents the partial derivative in the spatial variable $x$. The parameters ${\delta}, \widetilde{P}, P_e, D_a, K$ are positive real constants, while $\phi$ is a constant real number between 0 and 1. Their specific values, as they appear in ^[4], are summarized in Appendix A.1.

Following the standard derivations, we can obtain the following weak formulation for the above coupled system.

Definition 2.1. $(c, c_1, c_2)$ is called a weak solution for the system given by (2.1)–(2.10) if

and $c_2 \in H^1(0, T;L^2({\Omega^m}))$ satisfy that $c({\cdot}, 0) \equiv 1$ and $c_1({\cdot}, 0) \equiv c_2({\cdot}, 0) \equiv 0$, and, for some $T > 0$, any $t\in (0, T]$ and any $(v, w)\in H^1({\Omega^s})\times H^1({\Omega^m})$ such that

where $({\cdot}, {\cdot})$ denotes the $L^2$-inner product on the respective domain, $< {\cdot}, {\cdot} >$ denotes the duality pairing, and

with ${\delta}, \widetilde{P}, P_e, D_a, K$ being positive constant parameters and $\phi$ being a constant between 0 and 1.

For notation brevity, but without loss of clarity, throughout this paper, we may omit the explicit domain dependence in spatial norms. For example, $\| {\cdot}\|_{L^2({\Omega^s})}$ could be written as $\| {\cdot}\|_{L^2}$.

3.   PDE analysis

The goals of this section are to prove the well-posedness for the coupled PDE system given by (2.1)–(2.10) and establish some properties for the weak solution, including the boundedness property. To this end, we first derive some needed a priori estimates for weak solutions. Then, we prove the existence and uniqueness by using the Galerkin and energy methods.

3.1.   A priori estimates

The main result of this subsection is summarized in the following theorem.

Theorem 3.1. Let $(c, c_1, c_2)$ be a weak solution to the system given by (2.1)–(2.10) in the sense of Definition 2.1. Then, the following holds:

where $C > 0$ is a constant that is independent of $(c_1, c_2, c)$.

Proof. Setting $v: = c$ and $w: = c_1$ in Definition 2.1, by Theorem A.3, Hölder's inequality, and the fundamental theorem of calculus, we get

In order to cancel out boundary terms on the right-hand side, we decided to choose $\epsilon_1 = \epsilon_2 = \frac{1}{2}$ and $\epsilon_3 = \epsilon_4 = 1$ and let $\gamma: = \frac{1}{2} \min\{\phi, 1-\phi\}$. Combining (3.4) to (3.6) together, applying Lemma A.2 with

and taking the supremum over $[0, T]$ on both sides yields

which proves (3.1).

Next, to recover the boundary term $\|c(0, {\cdot})\|^2_{L^2(0, T)}$, notice that the above estimate, in particular, implies that $\|c_1(0, {\cdot})\|^2_{L^2(0, T)}$ is controlled from the above. Therefore, setting $\epsilon_1 = \frac{1}{2}$ and integrating (3.4) over $(0, T)$ yields

which gives (3.2).

We have yet to estimate the functional norms $\| {\partial}_t c\|_{H^{-1}({\Omega^s})}$ and $\| {\partial}_t c_1\|_{H^{-1}({\Omega^m})}$. It suffices to show that all terms in the bilinear forms are bounded. By Lemma A.1, we immediately have the following estimate:

Notice that the constant $C$ here depends on the spaces to which the functions $w$ and $v$ belong. If both $w, v\in H^1({\Omega^s})$, then $C = 2l^{-1}+1 = \mathcal{O}(l^{-1})$. If both $w, v\in H^1({\Omega^m})$, then $C = 3 = \mathcal{O}(1)$. If $w\in H^1({\Omega^s})$ and $v\in H^1({\Omega^m})$, or vice versa, then $C = (6l^{-1}+3)^{1/2} = \mathcal{O}(l^{-1/2})$, which is the only case that explicitly appears in Definition 2.1. However, the first two cases appear implicitly within the bilinear forms $\mathcal{A}[ {\cdot}, {\cdot}]$ and $\mathcal{B}[ {\cdot}, {\cdot}]$. Consequently, by Hölder's inequality, we get

Therefore,

where $C = C \big(l^{-1}, {\delta}, {\delta} \widetilde{P}, P_e, \phi^{-1}, (1-\phi)^{-1}, D_a, K^{-1} \big)$ with linear dependence on each argument. This then verifies (3.3) and hence concludes the proof. □

Remark 3.1. The boundedness now becomes a property of the weak solution via the compact embedding $H^1 ((x_1, x_2))\hookrightarrow L^\infty ((x_1, x_2))$ for any $x_1 < x_2$. This validates our modified model and the newly imposed no-flux boundary condition.

3.2.   Well-posedness

Since the system given by (2.1)–(2.10) is linear, uniqueness would be an immediate corollary of a priori estimates because, if there are two weak solutions, their difference must satisfy the conditions of the same equations but take zero initial conditions, which, in turn, implies that the difference is zero. Hence, we have the following theorem.

Theorem 3.2 (Uniqueness). There exists at most one weak solution $(c, c_1, c_2)$ to the problem given by (2.1)–(2.10) in the sense of Definition 2.1.

To prove the existence, we adopt the Galerkin method with a compactness argument. To setup our Galerkin approximation, let $\mathcal T^s : = \cup_{j = 1}^{N} I^s_j$ and $\mathcal T^m : = \cup_{j = 1}^{N} I^m_j$ be uniform meshes on ${\Omega^s}$ and ${\Omega^m}$ respectively. Let $\{\psi^s_j\}_{j = 1}^{N+1}, \{\psi^m_j\}_{j = 1}^{N+1}$ be the standard linear finite element nodal basis functions on $\mathcal T^s$ and $\mathcal T^m$, respectively, and define

where $h = l/N$ in $V^s_h$ and $h = 1/N$ in $V^m_h$. Here, we abuse the notation to give $h$ multiple meanings for the sake of notation brevity.

Definition 3.3. $(c_h, c_{1h}, c_{2h}): [0, T] \to V^s_h\times V^m_h\times V^m_h$ is called an approximate weak solution to the system given by (2.1)–(2.10) if the following holds for any $(v_h, w_h) \in V^s_h\times V^m_h$:

with the initial conditions $c_h({\cdot}, 0) \equiv 1$ and $c_{1h}({\cdot}, 0) \equiv c_{2h}({\cdot}, 0) \equiv 0$.

Lemma 3.4. For each $h > 0$, there exists a unique approximate weak solution $(c_h, c_{1h}, c_{2h})$ in the sense of Definition 3.3.

Proof. By the definition, $c_h$, $c_{1h}$, $c_{2h}$ can be written as follows:

Then, the equations in Definition 3.3 can be rewritten as follows:

for each $j = 1, \cdots, N+1$.

Equations (3.10) to (3.12) can be rewritten as the following ODE system:

where

and

We note that, for $j = 1, 2, \cdots, N+1$, the following holds:

Hence, the existence of a unique approximate weak solution is equivalent to the existence of a unique solution to the above ODE system. Since $\Psi_s$ and $\Psi_m$ are tri-diagonal and strictly diagonally dominant, they are invertible. Then, $\mathit{\boldsymbol{y}}\, ' (t) = {\mathit{\boldsymbol{D}}}^{-1} {\mathit{\boldsymbol{M}}} \, \mathit{\boldsymbol{y}}(t)$. Since ${\mathit{\boldsymbol{D}}}^{-1} {\mathit{\boldsymbol{M}}}$ is a constant matrix, by Lemma A.4, the ODE system has a unique solution. Thus, there exists a unique approximate weak solution $(c_h, c_{1h}, c_{2h})$. □

Theorem 3.5 (Existence). There exists a weak solution $(c, c_1, c_2)$ to the problem given by (2.1)–(2.10) in the sense of Definition 2.1.

Proof. We first notice that the approximate weak solution proved in Lemma 3.4 satisfies the conditions of those estimates of Theorem 3.1. Since $L^2(0, T;H^1({\Omega^s}))$ is a reflexive Banach space, $L^\infty(0, T;L^2({\Omega^s}))$ is a separable normed linear space, and the sequence $\{c_h\}$ is uniformly bounded in $h$, then there exists a subsequence $\{c_{h_j}\}$ that converges weakly and weak* in $L^\infty(0, T;L^2({\Omega^s}))$ and $L^\infty(0, T;L^2)$, respectively (see Theorems A.5 to A.7), that is, there exists $c\in L^2(0, T;H^1({\Omega^s})) \cap L^\infty(0, T;L^2({\Omega^s}))$ such that

Moreover, since $H^1(0, T;H^{-1}({\Omega^s}))$ is a separable normed linear space and $\{c'_h\}$ is uniformly bounded in $h$, there exists a subsequence of $\{c_{h_j}\}$ (not relabeling here for notation brevity) such that it converges in a weak* sense; namely, there exists $\zeta\in L^2(0, T;H^{-1}({\Omega^s}))$ such that

We want to show that $\zeta$ is actually the weak time derivative of $c$. To this end, let $\eta \in C_0^\infty(0, T)$ and $v_h \in V^s_h$; then,

Therefore, $\zeta$ is the weak time-derivative of $c$ by definition. Now, let $\eta \in C^\infty_0(0, T)$ and $v_h \in V^s_h$; by Definition 3.3, we have

For each fixed $\eta$ and $v_h$, notice that $\, \mathcal{A}[ {\cdot}, \eta v_h]$ defines a bounded linear functional on the space to which $\{c_{h_j}\}$ belongs. Therefore, setting $j\to \infty$, and by weak convergence, we get

Similarly,

Since $\eta\in C^\infty_0(0, T)$ is arbitrary, then the above equations infer that

which, with the denseness of $V^s_h$ in $H^1({\Omega^s})$, implies that

Hence, (2.11) holds.

Using exactly the same argument we can show that $c_1$ satisfies (2.12).

It remains to be shown that $c_2$ satisfies (2.13). To this end, let $\eta \in C_0^\infty (0, T)$ and $w_h\in V^m_h$; then, it follows that

Using a previous derivation, we can show that

In addition, since $\{c_{2h_j}\}$ is uniformly (in $h$) bounded in $L^\infty(0, T;L^2)$, which is a separable normed linear space, and $\{c'_{2h_j}\}$ is uniformly bounded in $L^2(0, T;L^2)$, which is a reflexive Banach space, then there exists a subsequence of $\{c_{2h_j}\}$ (not relabeling for notation brevity) such that, as $j\to\infty$,

Again, with the help of integration by parts and the definition, we can show that $\theta = {\partial}_t c_2$. Therefore, as $j\to \infty$,

Consequently,

Since $\eta\in C_0^\infty(0, T)$ and $w_h \in V_h^m$ are arbitrary and $V^m_h$ is dense in $H^1(\Omega^m)$, then

Thus, $(c, c_1, c_2)$ is a weak solution to the problem given by (2.1)–(2.10) by Definition 2.1. The proof is complete. □

Corollary 1 (Convergence). The finite-element approximate weak solution $(c_h, c_{1h}, c_{2h})$ converges to the unique PDE solution $(c, c_1, c_2)$.

Proof. From the proof of Theorem 3.5, we conclude that every convergent subsequence of the finite-element approximate weak solution $(c_h, c_{1h}, c_{2h})$ converges to a PDE weak solution $(c, c_1, c_2)$. Since the PDE weak solution is unique, the whole sequence $(c_h, c_{1h}, c_{2h})$ must converge to the PDE solution $(c, c_1, c_2)$. □

4.   Error estimates and formulation of fully discrete scheme

In the last section, we constructed a semi-discrete finite-element Galerkin approximation to the problem given by (2.1)–(2.10) and proved its convergence (see Corollary 1) as a byproduct of the proof of the existence theorem. The primary goals of this section are to derive optimal rates of convergence in powers of $h$ (i.e., error estimates) for the finite-element solution, and to formulate a practical fully discrete scheme which will be used in the subsequent section for numerical simulations and to numerically verify the sharpness of the proved convergence rates.

4.1.   Error estimates for semi-discrete finite element method

We recall that PDE and finite-element approximate solutions were respectively defined in Definitions 2.1 and 3.3. To derive the error estimates, we first need to obtain the error equations; to this end, subtracting the equations in Definition 3.3 from their corresponding equations in Definition 2.1 (with the same test functions $v_h \in V^s_h$ and $w_h \in V^m_h$), we get

where $e_h : = c-c_h, e_{1h} : = c_1-c_{1h}$, and $e_{2h} : = c_2-c_{2h}$.

Let $\mathcal{R}_h^A: H^1({\Omega^s}) \to V^s_h$ be the elliptic projection defined by

and $\mathcal{R}_h^B: H^1({\Omega^s}) \to V^m_h$ be another elliptic projection defined by

These projection operators are well-defined because each bilinear form is coercive and continuous. Further, let $\mathcal{P}_h$ be the $L^2$-projection onto $V^m_h$ defined by

We also introduce the following error decompositions:

It is well known (see ^[15]) that

Using the error decompositions we can rewrite the error equations (4.1)–(4.3) as follows:

To derive the desired error estimates, our task now is to control $\xi$ terms via the $\rho$ terms by using the above error equations.

Theorem 4.1 (Error estimates). The following error estimates hold:

Proof. Similar to the a-priori estimates, setting $\xi_h, \xi_{1h}$, and $\xi_{2h}$ to be the test functions, we obtain the following inequalities:

where

Choose appropriate values of $\epsilon_2$ and $\epsilon_5$ such that

Furthermore, choosing $\epsilon_1 = \epsilon_4 = \frac{1}{2}$, as well as

and adding (4.12) to (4.14), yields

where $E: = E_h + E_{1h} + E_{2h}$.

By using Gronwall's inequality given in Lemma A.2 and taking the supremum over $(0, T)$, we get

In particular,

By applying (4.12) with the above choices of $\epsilon_1$ and $\epsilon_2$, we have

Integrating over $(0, T)$ in $t$ yields

In particular,

In summary, we have shown that

which, combined with (4.4)–(4.6) and an application of the triangle inequality, concludes the proof. □

Remark 4.1. (a) Both estimates (4.10) and (4.11) are optimal compared to the linear finite-element interpolation errors.

(b) The $L^2$-norms of $(\xi_h)_x$ and $(\xi_{1h})_x$ exhibit a superconvergence property.

(c) If $r$th ($r > 1$)-order finite-element space is used in place of the linear finite-element space and we assume that the solution $(c, c_1, c_2)$ is sufficiently regular, then it can be proved that the rates of convergence in (4.10) and (4.11) will be improved to $O(h^{r+1})$ and $O(h^r)$, respectively.

4.2.   Formulation of fully discrete schemes

To get a computable fully discrete method, we need to discretize (3.10) to (3.12) in time by using any time-stepping scheme, such as the Euler, implicit Euler, Runge-Kutta, backward differentiation formula (BDF), and Crank-Nicolson methods. Below, we use the simplest Euler method to demonstrate the procedure. For each $j$, the Euler method is given by

They can be rewritten in matrix-vector form, as follows:

where ${\mathit{\boldsymbol{f}}}_0^k = [0;...; 0; y^k_{1, 1}]$ and ${\mathit{\boldsymbol{f}}}_1^k = [y^k_{0, N+1}; 0;...; 0]$.

It is well known that the Euler method results in an error of order $\mathcal{O}({\Delta} t)$; therefore, it can be shown that the fully discrete error is of order $\mathcal{O}({\Delta} t + h^2)$, provided that the Courant-Friedrichs-Lewy (CFL) condition ${\Delta} t < \min\{\frac{h^2}{2 {\delta}}, \frac{\phi h^2}{2}\}$ holds.

5.   Numerical simulations

In this section we present some numerical simulation results, which were computed by using Matlab R2022a. Since we aimed to solve a coupled system, a decoupling strategy was needed. In addition, due to the biological nature of the system, we propose a multi-rate time-stepping procedure to solve the system in order to save computation time. Several comparisons are made in Sections 5.1 to 5.2 among different schemes, different time-stepping strategies, and different decoupling strategies. Numerical results of the simulations are presented in Section 5.3.

The following notations are adopted in this section. Let $h$ denote the mesh size and ${\Delta} t$ the time step size. When those values differ in the two domains, we denote them as $h_ {\Omega^s}$, ${\Delta} t_ {\Omega^s}$ and $h_ {\Omega^m}$, ${\Delta} t_ {\Omega^m}$, respectively. Furthermore, let $N_0 : = l/h_ {\Omega^s}, N : = 1/h_ {\Omega^m}, n_0 : = T/ {\Delta} t_ {\Omega^s}$, and $n : = T/ {\Delta} t_ {\Omega^m}$.

5.1.   Comparison between two decoupling strategies

We propose two decoupling strategies for solving the system given by (4.15)–(4.17). The first strategy is parallelizable, updating $c$ and $c_2$ simultaneously, followed by updating $c_1$, as shown in Figure 3. The second strategy is a sequential update, which is shown in Figure 4.

Test cases were run for $T_{end} = 1, 10,100,200$, also applying different stepping strategies. In some cases, Algorithm I ran faster than Algorithm II, and, in some cases, it was vice versa. Furthermore, the difference in total CPU time was within 5% between the two algorithms. As for the accuracy comparison, we computed the numerical "true" solutions first, using $T = 1, \widetilde N_0 = 1000, \widetilde N = 500, \widetilde n = 2.5816\times 10^{6}.$ Then, the approximate solutions were computed for $T = 1, N_0 = 50, N = 25, n = 6454.$ The relative errors in $c$ and $c_2$ were both under $10^{-3}$; the relative error in $c_1$ was around $10^{-2}$ with Algorithm I, and it was about 10% more accurate than Algorithm II. This can be explained by the fact that $c_1$ is updated with the most recent coupling values in Algorithm I.

5.2.   Comparison between two time-stepping strategies

Since each subdomain possesses distinct biological/physical properties, it is natural to use different mesh sizes for different subdomains. In addition, the time step sizes can also be distinct in different subdomains; one scenario was that one time step size was taken as a constant multiple of the other. Recall that the Robin boundary condition of ${\Omega^s}$ at $x = 0$ states that $c_x + \widetilde{P} c = \widetilde{P} c_1$, where $\tilde{P} = 45000$. Therefore, it is natural to use a fine spatial mesh in ${\Omega^s}$. This does not have much effect on the mesh size despite the CFL condition since the diffusion coefficient was extremely small, at ${\delta} \sim 10^{-7}$. Therefore, the explicit time-stepping in fact demonstrates no disadvantage in this regard.

We restricted ourselves to Algorithm I and compared the performance of the naive and the multi-rate time-stepping strategies. When $N_0 = 2N$, all three errors decreased to around 10%, and, in the case of $c_2$, the errors decreased to around 1% of their naive counterparts, without increasing CPU time. Table 1 shows the relative errors of the three equations respectively.

It is worth pointing out that increasing the ratio $N_0/N$ further would increase CPU time without seeing any improvement in the accuracy. Furthermore, choosing $N_0 < N$ not only increased CPU time, it also resulted in larger error.

5.3.   Simulation results

We computed the concentrations $c, c_1, c_2$ with $T$ taken as 10 minutes, 30 minutes, 1 hour, 6 hours, and 24 hours after the stent insertion, using Algorithm I. Computed results are summarized in Figure 5 and Figure 6. It can be observed that the diffusion within the stent is slow. This is in line with the extremely small diffusion coefficient, i.e., ${\delta} \ll1$. However, at the interface, the advection of the drug into the media region is initially extremely fast since it is proportional to the concentration difference, which led to the steep drop near $x = 0$. As for the media concentrations, $c_1$ increased for a period of time before eventually dropping due to the absorption of the drug into the muscle cells, where a steady increase in drug concentration is demonstrated. In addition, Figure 7 shows the interface concentration $c_1(0, t)$ over the first 24 hours after stent insertion, reproducing the profile shape from ^[4], which was obtained via an analytic method.

The finite-element simulation results presented in this section were also confirmed by those obtained via a finite-difference method (which are not presented here to save space). The $L^\infty(L^2)$-error between the finite-element and the finite-difference results was found to be around $5 \times 10^{-5}$ when $h_ {\Omega^s} = 2.8 \times 10^{-4}$ and $h_ {\Omega^m} = 0.01$.

6.   Two-dimensional generalization

Below, we propose a two-dimensional drug-release model which is an extension of the one-dimensional model studied in the previous sections. Figure 8 shows a two-dimensional sketch of the arterial wall.

Again, as in the one-dimensional case, $c, c_1,$ and $c_2$ represent the unknown concentrations in the stent, in the extracellular matrix, and in the muscle cells, respectively. All other quantities are given coefficients. Unlike the one-dimensional model, we note that there are four additional pieces of the boundary. New boundary conditions must be prescribed there. For these four pieces, Dirichlet boundary conditions were imposed on $\Gamma^a$ and $\Gamma^b$, whereas the periodic boundary conditions were used on $\Gamma^c$ and $\Gamma^d$ for $c_1$.

Essentially, all of the PDE and numerical analysis results for the 1-d model can be extended to the 2-d model; this includes a priori estimates and well-posedness proofs, as well as the finite-element convergence and error estimates. One notable difference is that the boundary of the 2-d domain consists of 1-d line segments; handling the 2-d boundary terms requires some additional and more involved trace inequalities. Below, we only state the key formulations and main results, without giving detailed derivations and proofs, to save space. However, it should be noted that, although the PDE and numerical analyses are similar, the computer simulations and coding are much more complicated in the 2-d case; we will present those results elsewhere.

Our two-dimensional model is given as follows:

We note that (6.10) and (6.11) represent the periodic boundary conditions for $c_1$.

Theorem 6.1 (Well-posedness in 2-d). Under some reasonable assumptions on the coefficients, there exists a unique weak solution $(c, c_1, c_2)$ satisfying the following for any $T > 0$:

and, for $t\in (0, T]$,

Here,

$({\cdot}, {\cdot})$ denotes the $L^2$-inner product, $< {\cdot}, {\cdot} >$ denotes the duality pairing, and

Before introducing the finite-element result, we first introduce the finite-element spaces and the concept of approximate weak solutions. Let $\mathcal T^s_h$ be a mesh of ${\Omega^s}$ and $\mathcal T^m_h$ a mesh of ${\Omega^m}$, both geometrically conformal. Let

It is a known fact that each $\mathbb V_h^s$ and $\mathbb V_h^m$ has a set of nodal basis functions, denoted by $\{\psi_j^s\}$ and $\{\psi_j^m\}$, respectively, and satisfying that $\psi_j^s({\mathit{\boldsymbol{a}}}_i) = {\delta}_{ij}$ and $\psi_j^m(\hat{{\mathit{\boldsymbol{a}}}_i}) = {\delta}_{ij}$ for each node ${\mathit{\boldsymbol{a}}}_i$ and $\hat{{\mathit{\boldsymbol{a}}}_i}$, and each $j$.

Definition 6.2. $(c_h, c_{1h}, c_{2h}): (0, T] \rightarrow \mathbb V^s_h \times \mathbb V^m_{h, per} \times \mathbb V^m_h$ is called an approximate weak solution to the 2-d system if

with $c_h(0, {\cdot}) = \mathcal{P}_h c_0 \in \mathbb V^s_h$ and $c_{1h}(0, {\cdot}) \equiv c_{2h}(0, {\cdot}) \equiv 0$.

We obtained the following well-posedness and error estimate results.

Theorem 6.3 (Error estimates in 2-d). For each $h > 0$, there exists a unique approximate weak solution $(c_h, c_{1h}, c_{2h})$. Moreover, let $(c, c_1, c_2)$ be the weak solution stated in Theorem 6.1 and define error functions $e_h : = c-c_h$, $e_{1h} : = c_1-c_{1h}$, and $e_{2h} : = c_2-c_{2h}$, then, the following error estimates hold:

We note that the interface terms now appear as $L^2$ integrals on $\Gamma$ in the 2-d case, which only gives us the control of $\|e_h\|_{L^2(0, T;L^2(\Gamma))}$ and $\|e_{1h}\|_{L^2(0, T; L^2(\Gamma))}$, not pointwise control as in the 1-d case. But, these estimates are consequences of (6.16) and the trace inequality.

7.   Summary and concluding remarks

In this paper we have presented a complete PDE and numerical analysis for the modified one-dimensional intravascular stent model originally proposed in ^[4]. The modified model is not only well-posed, it also has improved numerical computability. The well-posedness was proved by using the Galerkin method in combination with a compactness argument. A semi-discrete finite-element method and a fully discrete scheme using the Euler time-stepping was formulated for the PDE model. Optimal order error estimates in the energy norm was proved for both schemes. Numerical results have been presented, along with comparisons between different decoupling strategies and time-stepping schemes. Finally, extensions of the model, along with its PDE and numerical analysis results for the two-dimensional case, were also briefly discussed.

Our future work on this research project will focus on the direct generalizations of this model in one and more dimensions in the modeling of the transmural advection using Darcy's flow; the model will also include an additional lumenal layer that considers the effect of blood flow on drug delivery. Moreover, we are very much interested in completing higher-dimensional codes and simulations as well as in validating the model simulations with real lab data.

It is also worth noting that the analysis techniques employed in this work should be readily extendable to other more sophisticated linear models. In addition, there have been recent developments in nonlinear drug-binding models (see ^[16] and ^[17]). In those cases, the nonlinear terms require special care and new techniques. It is expected that our general techniques might still be applicable. Such a detailed analysis will be carried out as future work.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

This work was partially supported by the NSF grants DMS-2012414 and DMS-2309626.

Conflict of interest

The authors declare that there is no conflict of interest.

A.   Parameters and cited facts

A.1.   Parameters

We summarize in the chart below the parameter values that appeared in the one-dimensional model and used in our simulations. We refer the reader to ^[4] for more details.

A.2.   Some known facts

In this subsection, we present a few well-known lemmas and theorems cited in this paper and either provide a brief proof or cite at least one reference where a proof can be found.

Lemma A.1 (A 1-d trace inequality). Let $v: [a, b]\to \mathbb{R}$, $v\in H^1(a, b)$, and $H: = b-a$. Then

where $C \approx 1$ if $H\gg1$ and $C \in \mathcal{O} \big(H^{-1/2} \big)$ if $H\ll1$.

Proof. Let $w: [a, b]\to \mathbb{R}$ be differentiable. By the fundamental theorem of calculus we have the following for any $a\leq w_1 < w_2\leq b$:

Setting $w(x): = [(x-a)v(x)]^2$, $x_1: = a$, and $x_2: = b$ yields

Therefore,

which gives the first inequality. The second inequality follows similarly with $w(x): = [(b-x)v(x)]^2$. □

Lemma A.2 (General form of Grönwall's inequality; see Appendix B2 of ^[18]). Let $\xi, \phi, \psi \in L^2(0, T)$ and be nonnegative, and let $\eta \in H^1(0, T)$. If

then

Theorem A.3 (Theorem 5.9 of ^[18]). Suppose that ${u} \in L^2(0, T;H^1_0(U))$ and ${u}' \in L^2(0, T;H^{-1}(U))$.

(i) Then, ${u} \in C([0, T];L^2(U))$ (after possibly being redefined on a set of measure zero).

(ii) The mapping $t \mapsto \|{u}(t)\|^2_{L^2(U)}$ is absolutely continuous, with

(iii) Furthermore, the following holds

Lemma A.4 (Lemma 1.1, Chapter IV of ^[19]). The initial value problem given by

has a unique solution $\mathit{\boldsymbol{y}}:[0, T]\to \mathit{\boldsymbol{R}}^n$ if $Q$ is a continuous function on $[0, T]$.

Theorem A.5 (Theorem 4.10.8 of ^[20]). Let $X$ be a reflexive Banach space. A set $K\subset X$ is weakly sequentially compact if and only if it is both bounded and weakly closed.

Theorem A.6 (Theorem 4.12.3 of ^[20]). Let $X$ be a separable normed linear space. Then every bounded sequence of continuous linear functionals in $X^*$ has a weakly convergent subsequence.

Theorem A.7 (Theorem 5.1.1 of ^[20]). A completely continuous (compact) linear operator maps weakly convergent sequences into convergent sequences.