# Quadratic ODE and PDE Models of Drug Release Kinetics from Biodegradable Polymers

Michel C. Delfour, André Garon

► To cite this version:

Michel C. Delfour, André Garon. Quadratic ODE and PDE Models of Drug Release Kinetics from Biodegradable Polymers. 25th System Modeling and Optimization (CSMO), Sep 2011, Berlin, Germany. pp.13-24, 10.1007/978-3-642-36062-6_2. hal-01347516

HAL Id: hal-01347516

https://inria.hal.science/hal-01347516

Submitted on 21 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Quadratic ODE and PDE models of drug release kinetics from biodegradable polymers

Michel C. Delfour\* and André Garon \*\*

delfour@crm.umontreal.ca
Andre.Garon@polymtl.ca

**Abstract.** In order to achieve prescribed drug release kinetics over long therapeutic periods, bi-phasic and possibly multi-phasic releases from blends of biodegradable polymers are currently envisioned. The modelling of drug release in the presence of degradation of the polymer matrix and surface erosion is quite complex. Yet, simple reliable mathematical models validated against experimental data are now available to help in classifying neat polymers and in predicting the release dynamics from polymer blends. In this paper, we survey a two-parameter quadratic ODE model that has been validated against experimental data for the release of paclitaxel from a broad range of biodegradable polymers and a quadratic semi-permeable membrane PDE model that mimics the ODE model and could readily be extended to drug eluting stents.

**Keywords:** Drug release models, biodegradable polymers, paclitaxel.

# 1 Introduction

Stents are used in interventional cardiology to keep a diseased vessel open after angioplasty. This procedure is known to damage the endothelium at the insertion site and thus to favour the occurrence of in-stent restenosis through the proliferation of smooth muscle cells (SMC) within the vessel lumen. To control the abnormal behaviour of SMC, stents are coated with polymers that slowly release drug through diffusion into the vessel wall (drug-eluting stents or DES). These drugs are designed to control the rate of mitosis of SMC until the regeneration of the endothelium. The reader is referred to T. Kataoka et als [15] in 2002 and Joner et als [14] in 2006 for a fairly well-documented account of DES for the prevention of neointimal growth (see, for instance, [15, Figure 1, p. 1791]).

If endothelial cells do not recover to effectively control the proliferation of SMC's, a sustained dose will be required over the therapeutic period and even forever. In order to achieve prescribed drug release kinetics the current design

<footnote>\* Centre de recherches mathématiques et Département de mathématiques et de statistique, Université de Montréal, CP 6128, succ Centre-ville, Montréal (Qc), Canada H3C 3J7.</footnote>

<footnote>\*\* Département de Génie mécanique, École Polytechnique de Montréal, C.P. 6079, succ. Centre-ville, Montréal (Qc), Canada H3C 3A7.</footnote>
strategies focus on bi-phasic and possibly multi-phasic releases from blends of
biodegradable polymers (see, for instance, Batycky et al [1] in 1997) to achieve
specific drug release kinetics profiles over long therapeutic windows.

Recently, Lao and Venkatraman [16] and Lao, Venkatraman, and Peppas [18] have proposed a semi-empirical model to predict the release profile of paclitaxel from three neat polymer matrices: PCL (Polycaprolactone), PLGA (dl-lactide-co-glycolide) and PLGAPEG (PLGA with polyethylene glycol). They are representative of a broad family of biodegradable polymers ranging from hydrophobic to hydrophilic. In hydrophilic polymers the internal bounds between the chains are weakened and this adds to the surface erosion phenomenon. The drug release mechanism within a polymer matrix depends on many factors such as the affinity of the drug with the surrounding medium (water). Specifically, paclitaxel is hydrophobic and this might explain the fact that some of the drug blended into the polymer matrix is not released and cannot participate to the treatment of the disease wall. This is a difficult subject. The main criticism expressed in [18] of available models for drug release from eroding surfaces is that they fail to faithfully reproduce experimental data for highly degradable polymers (the S-curve behaviour). The reader is referred to the introduction of the paper of Lao et al [18] for a comprehensive review of the literature.

A quick look at the paclitaxel release profiles suggests two types of release: S-curve type and exponential type. S-curve behaviours are similar to the ones encountered in the study of the logistic equation of populations. In [2] we introduced a simple two-parameter Ordinary Differential Equation (ODE) model that completely describes the paclitaxel release profiles from neat PCL, PLGA PEG, and PLGA polymer matrices. This model describes with greater accuracy the drug-release than the semi-empirical model of Lao et als [18] using 5 to 8 parameters.

The simplicity of our model for such a broad range of polymers indicates
that somehow the quadratic structure captures the complex microphysics and
chemistry of the release and degradation processes. Using a purely mathematical
intuition to modelling, we have introduced in [6] a time-space three dimensional
partial differential equation (PDE) model of the paclitaxel release that mimics
the ODE model. The film of neat polymer is modelled as a thin flat domain
whose polymer/medium interface is a quadratic semi-permeable membrane with
a concentration jump at the interface.

In this approach, the diffusion process through a semi-permeable membrane is modelled as a diffusion through an interface with cracks (not to be confused with holes) where the rate of transfer of the product is proportional to the size of the concentration jump across the interface. Since the cracks have zero surface, their *size* is measured in terms of the mathematical notion of *capacity*. What is very nice about this approach is that it is based on a mathematically well documented linear model coming from the study of the Neumann sieve by Damlamian [5] in 1985. It provides a variational formulation and a mathematically tractable approach to the asymptotic analysis of a punctured membrane as the size of the holes goes to zero while preserving a strictly positive capacity that accounts for
the diffusion of the drug through the cracks<fnref n="1" /> Adding the non-linearity captures the effect of the internal degradation of the polymer by making the rate of mass transfer proportional to the size of the concentration jump across the interface.

Our approach is different from others in the literature since it deals with the nonlinearity through a quadratic condition at the interface between the polymer and the surrounding medium instead of using a time-dependent or a nonlinear diffusion. This model can be seen as a first step towards a three dimensional modelling of the release of paclitaxel from drug eluding stents coated with biodegradable polymers. It is capable of covering a wide range of biodegradable polymers potentially including the ones for which an incomplete release is experimentally observed (recall that the paclitaxel is hydrophobic).

To complete the experimental approach to this modelling, the next step would be to set up an experimental benchmark to check if the model and the mathematical assumptions on the coefficients of the model are realistic. The validation of such a model would improve the modelling of the drug release part of the global three-dimensional model of a blood vessel incorporating the lumen, the blood, the aggregated wall, and the coated stent (cf. for instance, [8]) and the subsequent studies of the effect of the pattern of the stent in [3] and the effect of the pulsative nature of the blood in [7]. Such global studies are important to determine the set of features in the modelling of the blood vessel and of the stent that should be retained in the design of the stent and the drug release dynamics.

## 2 ODE model and gradient flow interpretation

In the previous paper [2] we have shown an excellent fit between experimental release data [16] of paclitaxel from biodegradable neat polymers and a two-parameter quadratic ODE model of the Riccati type. We briefly recall this model.

Given an initial mass $M_0 > 0$ of drug uniformly impregnated into a *polymeric matrix*, denote by $M(t) > 0$ the *mass of drug released outside the polymer* as a function of the time $t > 0$. Denote by $M_\infty, 0 \le M_\infty \le M_0$, the asymptotic mass of the drug released. The ODE model was chosen of the form

$$ \frac{dM}{dt}(t) = h(M(t)), \quad t > 0, \quad M(0) = 0, \qquad (2.1) $$

for some quadratic right-hand side

$$ h(M) \stackrel{\text{def}}{=} A_1 (M_\infty - M(t)) + A_2 (M_\infty - M(t))^2 \qquad (2.2) $$

such that $M'(0) = (A_1 + A_2 M_\infty) M_\infty > 0$. By introducing the *normalized released mass*

$$ m(t) \stackrel{\text{def}}{=} M(t)/M_0, \qquad (2.3) $$

<footnote n="1">¹ See also the more recent comprehensive paper [4, Theorem 5.5] using the very nice theory of periodic unfolding.</footnote>
we get the following quadratic ODE model

$$ \frac{dm}{dt}(t) = \left[ A_1 + A_2 M_0 \left( \frac{M_\infty}{M_0} - m(t) \right) \right] \left( \frac{M_\infty}{M_0} - m(t) \right), \quad m(0) = 0. \quad (2.4) $$

Assuming that the ratio $0 < M_\infty/M_0 \le 1$ is known, the model is completely specified by the two parameters $A_1$ and $A_2 M_\infty$. When $A_2 = 0$, the model is linear; when $A_2 \ne 0$, the right-hand side is of the form

$$ h(m) \stackrel{\text{def}}{=} A_2 M_0 (m_2 - m)(m_1 - m), \quad m_1 \stackrel{\text{def}}{=} \frac{M_\infty}{M_0}, \quad m_2 \stackrel{\text{def}}{=} \frac{A_1 + A_2 M_\infty}{A_2 M_0}. $$

It was shown in [2] that the following four cases can occur under the conditions $m(0) = 0$ and $m'(0) = A_1 + A_2 M_\infty > 0$:

Case 1) (True S type)

$$ A_1 > 0, A_2 < 0, \text{ and } -m_1 < \frac{1}{2} \frac{A_1}{A_2 M_0} \quad (\text{that is, } -m_1 < m_2 < 0), \quad (2.5) $$

with solution

$$ m(t) = m_1 m_2 \frac{1 - e^{-A_1 t}}{m_2 - m_1 e^{-A_1 t}} $$

for which the point of inflexion occurs at time $t_c = -(\log(-m_2/m_1))/A_1 > 0$;

Case 2) (S type)

$$ A_1 > 0, A_2 < 0, \text{ and } \frac{1}{2} \frac{A_1}{A_2 M_0} \le -m_1 \quad (\text{that is, } m_2 \le -m_1), \quad (2.6) $$

with the solution and the point of inflexion

$$ m(t) = m_1 m_2 \frac{1 - e^{-A_1 t}}{m_2 - m_1 e^{-A_1 t}}, \quad t_c = -(\log(-m_2/m_1))/A_1 \le 0; $$

Case 3) (Exponential type)

$$ A_1 \ge 0 \text{ and } A_2 > 0 \quad (\text{that is, } m_2 \ge 1), \quad (2.7) $$

with the solution and the blow up time

$$ \left\{ \begin{array}{l} m(t) = m_1 m_2 \frac{1 - e^{-A_1 t}}{m_2 - m_1 e^{-A_1 t}}, \quad \text{for } A_1 > 0 \text{ since } m_2 > m_1, \\ t_c = -\frac{\log(m_2/m_1)}{A_1} < 0, \end{array} \right. \quad (2.8) $$
$$ \left\{ \begin{array}{l} m(t) = m_1 \frac{A_2 M_\infty t}{1 + A_2 M_\infty t}, \quad \text{for } A_1 = 0 \text{ since } m_2 = m_1; \\ t_c = -\frac{1}{A_2 M_\infty} < 0, \end{array} \right. $$
Case 4) (True exponential) $A_1 > 0$ and $A_2 = 0$ with the solution

$$m(t) = \frac{M_{\infty}}{M_0} (1 - e^{-A_1 t}), \quad t_c = -\infty.$$

The generic behaviours of the solution $m(t)$ in the above cases are illustrated in Figures 1 and 2 and the parameters tabulated in Table 1 of [6].

The ODE model has an interesting gradient flow interpretation by introducing the function

$$E(m) \stackrel{\text{def}}{=} \frac{1}{2} A_1 \left( \frac{M_{\infty}}{M_0} - m \right)^2 + \frac{1}{3} A_2 M_0 \left( \frac{M_{\infty}}{M_0} - m \right)^3 \quad (2.9)$$

with gradient (derivative)

$$E'(m) = -A_1 \left( \frac{M_{\infty}}{M_0} - m \right) - A_2 M_0 \left( \frac{M_{\infty}}{M_0} - m \right)^2 \quad (2.10)$$

and Hessian (second order derivative)

$$E''(m) = A_1 + 2 A_2 M_0 \left( \frac{M_{\infty}}{M_0} - m \right). \quad (2.11)$$

The ODE can now be rewritten in the form of a gradient flow equation

$$\frac{dm}{dt}(t) + E'(m(t)) = 0, \quad m(0) = 0. \quad (2.12)$$

This is the continuous version of a steepest descent method to minimize the functional $E$. So, it is expected that starting from $m(0) = 0$ with $m'(0) > 0$ the asymptotic value $m_1$ of the solution of the ODE (2.12) would achieve a local minimum of $E(m)$. To do that, we compute the second derivative of $E$ under the assumption that $A_1 \ge 0$ and $m'(0) > 0$ which is equivalent to $E'(0) = -(M_{\infty}/M_0)[A_1 + A_2 M_{\infty}] < 0$. It turns out that in all cases except the second part of case 3), $m_1$ is a local minimum of $E(m)$. The exception corresponds to a point of inflection that can be changed into a global minimum by modifying the function $E$ to $E(m) = (A_2 M_0/3) |M_{\infty}/M_0 - m|^3$.

# 3 PDE model of quadratic semi-permeable membranes

## 3.1 Equations in the polymer and the surrounding medium

The experimental benchmark of [16] is contained in a vial. The polymer film is deposited flat at the bottom of the vial and the vial is filled with a fluid that we shall call the *surrounding medium* (see Figure 1). The vial is closed without circulation of the fluid. Denote by $\Omega_p$ the open domain occupied by the polymer and by $\Omega_m$ the open domain occupied by the surrounding medium. Let $\Gamma_p$ and $\Gamma_m$ be the respective boundaries of $\Omega_p$ and $\Omega_m$. The polymer occupies a thin square parallelepipedic region at the bottom of the vial. Its boundary is made up of the *interface* $\Gamma_{int} = \Gamma_p \cap \Gamma_m$ between the polymer and the medium (top
<figure><img src="image_2.png" /><figcaption>Fig. 1. The polymer film and the surrounding fluid in the vial.</figcaption></figure>

boundary and lateral boundary of $\Omega_p$) and the bottom square boundary of $\Omega_p$ that we shall denote $\Gamma_0$. The vial is closed without circulation of the fluid filling the vial (the medium). Within the (surrounding) medium only linear diffusion is expected with zero Neumann boundary conditions at the boundary of the vial $\Gamma_{ext} = (\Gamma_p \cup \Gamma_m)\setminus\Gamma_{int}$.

At time $t$, denote by $c_p(x, t)$ the concentration of the drug at the point $x \in \Omega_p$ and by $c_m(x, t)$ the concentration of the drug at point $x \in \Omega_m$. Assume linear diffusion equations in the polymer and the surrounding medium

$$ \frac{\partial c_p}{\partial t} = \text{div}(D_p \nabla c_p) \text{ in } \Omega_p \quad (3.1) $$

$$ \frac{\partial c_m}{\partial t} = \text{div}(D_m \nabla c_m) \text{ in } \Omega_m \quad (3.2) $$

with constant diffusion constants $D_p$ and $D_m$ and initial conditions

$$ c_p(x, 0) = c_0(x) = M_0/|\Omega_p| \text{ in } \Omega_p, \quad c_m(x, 0) = 0 \text{ in } \Omega_m, \quad (3.3) $$

where $|\Omega_p|$ is the volume of $\Omega_p$. Assume that the experimental set up is closed:

$$ D_p \frac{\partial c_p}{\partial n_p} = 0 \text{eq.constrictant} \Gamma_p \setminus \Gamma_{int} \quad D_m \frac{\partial c_m}{\partial n_m} = 0 \text{ on } \Gamma_m \setminus \Gamma_{int}, \quad (3.4) $$

where the unit normals $n_p$ and $n_m$ are exterior to the respective domains $\Omega_p$ and $\Omega_m$. Assume that there is no loss of product: this yields the (affine) constraint

$$ \forall t \ge 0, \quad M_0 \stackrel{\text{def}}{=} \int_{\Omega_p} c_0(x) \, dx = \int_{\Omega_p} c_p(x, t) \, dx + \int_{\Omega_m} c_m(x, t) \, dx, \quad (3.5) $$

where $M_0$ is the total mass of product. By integrating (3.1) over $\Omega_p$ and (3.2) over $\Omega_m$ and by using the constraint (3.5), we get

$$ \begin{aligned} \Rightarrow 0 &= \int_{\Omega_p} \frac{\partial c_p}{\partial t}(x, t) \, dx + \int_{\Omega_m} \frac{\partial c_m}{\partial t}(x, t) \, dx \\ &= \int_{\Omega_p} \text{div}(D_p \nabla c_p)(x, t) \, dx + \int_{\Omega_m} \text{div}(D_m \nabla c_m)(x, t) \, dx \quad (3.6) \\ &= \int_{\Gamma_p} D_p \frac{\partial c_p}{\partial n_p}(x, t) \, d\Gamma + \int_{\Gamma_m} D_m \frac{\partial c_m}{\partial n_m}(x, t) \, d\Gamma. \end{aligned} $$
Finally, by using the boundary conditions (3.4) we get

$$ \int_{\Gamma_{int}} \left[ D_p \frac{\partial c_p}{\partial n_p}(x,t) + D_m \frac{\partial c_m}{\partial n_m}(x,t) \right] d\Gamma = 0. \quad (3.7) $$

It remains to specify the conditions at the interface $\Gamma_{int}$.

## 3.2 Conditions at the interface

In order to incorporate the microphysics taking place in the thin film of polymer, it is assumed that the interface behaves as a semi-permeable membrane with micro fissures through which the drug diffuses into the surrounding medium. Many empirical and theoretical models of such membranes have been studied in the literature and in different contexts. One mathematically interesting model of a semi-permeable membrane is to assume that the interface is a membrane punctured with small holes whose size goes to zero while preserving a strictly positive capacity<fnref n="2" /> in the limiting process. In other words the membrane is fissured or cracked and the drug diffuses through the cracks. This problem has been studied from the mathematical point of view under the name of the Neumann sieve by A. Damlamian [5] in 1985. From the physical point of view, it can be assimilated with a semi-permeable membrane.

In this section we consider an evolution equation of the form

$$ \frac{\partial c}{\partial t}(t) + A(c(t)) = 0, \quad c(0) = M_0/|\Omega_p| \chi_{\Omega_p}, \quad (3.8) $$

where the operator A is now quadratic in $c(t)$. Since the domain $\Omega_p$ is thin, it is reasonable to put the nonlinearity at the interface $\Gamma_{int}$ rather than on $\Omega_p$ via a diffusion coefficient $D_p(c)$ that depends on c:

$$ \begin{aligned} -\frac{d}{dt} \int_{\Omega_p} c_p(t) dx &= \frac{d}{dt} \int_{\Omega_m} c_m(t) dx \\ &= \int_{\Gamma_{int}} \left[ k_1 + k_2 \frac{|\Omega_p|}{M_0} |c_p(t) - c_m(t)| \right] (c_p(t) - c_m(t)) d\Gamma \end{aligned} \quad (3.9) $$

for some constant $k_2$. Note that we have introduced a scaling by the initial concentration of product $M_0/|\Omega_p|$ of the drug so that $k_1$ and $k_2$ are parameters of the same physical dimension.

Now consider the (cubic) functional

$$ \begin{aligned} E(v) &\stackrel{\text{def}}{=} \frac{1}{2} \int_{\Omega_p} D_p |\nabla v_p|^2 dx + \frac{1}{2} \int_{\Omega_m} D_m |\nabla v_m|^2 dx \\ &\quad + \int_{\Gamma_{int}} \left( \frac{1}{2} k_1 |v_p - v_m|^2 + \frac{1}{3} k_2 \frac{|\Omega_p|}{M_0} |v_p - v_m|^3 \right) d\Gamma \end{aligned} \quad (3.10) $$

$$ v_p \stackrel{\text{def}}{=} v|_{\Omega_p}, \quad v_m \stackrel{\text{def}}{=} v|_{\Omega_m} \quad (3.11) $$

<footnote n="2">² The capacity of a set is a mathematical notion. For instance a finite segment in the plane has zero area but finite capacity. Roughly speaking, the capacity is a “measure” of the cracks.</footnote>
defined on the space $H^1(\Omega_p \cup \Omega_m)$ with a crack $\Gamma_{int} = \Gamma_p \cap \Gamma_m$ in $\Omega_p \cup \Omega_m$ along which the function $v$ can have a jump discontinuity $[v] = v_m - v_p$. This convex non quadratic variational formulation is similar to the $T^4$ radiation law for the temperature $T$ of a radiating body in free space (cf., for instance, [9]).

We do not impose the continuity of the concentrations at the interface. Taking into account the constraint on the total mass of product, we look for a solution $c(t)$ at time $t > 0$ in the affine subspace

$$V_{M_0}^{pm} \stackrel{\text{def}}{=} \left\{ \begin{array}{ll} \left\{ c \in H^1(\Omega_p \cup \Omega_m) : \int_{\Omega_p \cup \Omega_m} c(x) \, dx = M_0 \right\}, & \text{if } k_2 = 0, k_1 > 0, \\ \\ \left\{ c \in H^1(\Omega_p \cup \Omega_m) : \begin{array}{l} \int_{\Omega_p \cup \Omega_m} c(x) \, dx = M_0 \\ v_p - v_m \in L^3(\Gamma_{int}) \end{array} \right\}, & \text{if } k_2 > 0, \end{array} \right.$$

of $H^1(\Omega_p \cup \Omega_m)$. In the first case

$$\left[ \int_{\Omega_p} |\nabla v_p|^2 dx + \int_{\Omega_m} |\nabla v_m|^2 dx + \int_{\Gamma_{int}} |v_p - v_m|^2 d\Gamma \right]^{1/2} \quad (3.12)$$

is an equivalent norm on $V_{M_0}^{pm}$; in the second case

$$\left[ \int_{\Omega_p} |\nabla v_p|^2 dx + \int_{\Omega_m} |\nabla v_m|^2 dx \right]^{1/2} + \left[ \int_{\Gamma_{int}} |v_p - v_m|^3 d\Gamma \right]^{1/3} \quad (3.13)$$

is an equivalent norm on $V_{M_0}^{pm}$ (cf., for instance, [9] with the $T^4$ radiation law for the temperature $T$ in free space).

The directional derivative of $E$ is

$$\begin{aligned} dE(u; v) = & \int_{\Omega_p} D_p \nabla u \cdot \nabla v \, dx + \int_{\Omega_m} D_m \nabla u \cdot \nabla v \, dx \\ & + \int_{\Gamma_{int}} k_2 \frac{|\Omega_p|}{M_0} |u_p - u_m| (u_p - u_m) (v_p - v_m) \\ & + k_1 (u_p - u_m) (v_p - v_m) d\Gamma \end{aligned} \quad (3.14)$$

$$u_p \stackrel{\text{def}}{=} u|_{\Omega_p}, \quad u_m \stackrel{\text{def}}{=} u|_{\Omega_m}, \quad v_p \stackrel{\text{def}}{=} v|_{\Omega_p}, \quad v_m \stackrel{\text{def}}{=} v|_{\Omega_m}. \quad (3.15)$$

We are interested in the stationary points $c = (c_p, c_m) \in V_{M_0}$ of $E$ that are the solutions of the variational equation

$$\exists c \in V_{M_0}^{pm}, \quad dE(c; v) = 0, \quad \forall v \in V_0^{pm}, \quad (3.16)$$

$$V_0^{pm} \stackrel{\text{def}}{=} \left\{ \begin{array}{ll} \left\{ c \in H^1(\Omega_p \cup \Omega_m) : \int_{\Omega_p \cup \Omega_m} c(x) \, dx = 0 \right\}, & \text{if } k_2 = 0, k_1 > 0, \\ \\ \left\{ c \in H^1(\Omega_p \cup \Omega_m) : \begin{array}{l} \int_{\Omega_p \cup \Omega_m} c(x) \, dx = 0 \\ v_p - v_m \in L^3(\Gamma_{int}) \end{array} \right\}, & \text{if } k_2 > 0. \end{array} \right. \quad (3.17)$$
Again $dE(c; v) = 0$ for all constant functions $v$ and $V_0^{pm}$ can be replaced by $H^1(\Omega_p \cup \Omega_m)$:

$$ \exists c \in V_{M_0}^{pm}, \quad dE(c; v) = 0, \quad \forall v \in H^1(\Omega_p \cup \Omega_m). $$

It yields a complete set of conditions at the interface and the following system of equations

$$ \operatorname{div} (D_p \nabla c_p) = 0 \text{ in } \Omega_p, \quad \operatorname{div} (D_p \nabla c_m) = 0 \text{ in } \Omega_m \quad (3.18) $$

$$ D_p \frac{\partial c_p}{\partial n_p} + k_2 \frac{|\Omega_p|}{M_0} |c_p - c_m| (c_p - c_m) + k_1 (c_p - c_m) = 0 \text{ on } \Gamma_{int} \quad (3.19) $$

$$ D_m \frac{\partial c_m}{\partial n_m} - \left[ k_2 \frac{|\Omega_p|}{M_0} |c_p - c_m| (c_p - c_m) + k_1 (c_p - c_m) \right] = 0 \text{ on } \Gamma_{int} \quad (3.20) $$

$$ D_p \frac{\partial c_p}{\partial n_p} = 0 \text{ on } \Gamma_p \setminus \Gamma_{int}, \quad D_m \frac{\partial c_m}{\partial n_m} = 0 \text{ on } \Gamma_m \setminus \Gamma_{int} \quad (3.21) $$

$$ \int_{\Omega_p} c_p \, dx + \int_{\Omega_m} c_m \, dx = M_0. \quad (3.22) $$

From the mathematical viewpoint, the condition involving $|c_p - c_m| (c_p - c_m)$ is the analogue of the condition $|T - T_m|^3 (T - T_m)$ (usually written $(T - T_m)^4$) on the temperature of a radiating body (cf., for instance, [9]). The thin layer of polymer behaves as a nonlinear semi-permeable membrane. The second order directional derivative of $E$ is

$$ d^2 E(u; v; w) = \int_{\Omega_p} D_p \nabla w \cdot \nabla v \, dx + \int_{\Omega_m} D_m \nabla w \cdot \nabla v \, dx \\ + \int_{\Gamma_{int}} \left[ 2 k_2 \frac{|\Omega_p|}{M_0} |u_p - u_m| + k_1 \right] (w_p - w_m) (v_p - v_m) \, d\Gamma \quad (3.23) $$

$$ \Rightarrow \ d^2 E(u; v; v) = \int_{\Omega_p} D_p |\nabla v|^2 \, dx + \int_{\Omega_m} D_m |\nabla v|^2 \, dx \\ + \int_{\Gamma_{int}} \left[ 2 k_2 \frac{|\Omega_p|}{M_0} |u_p - u_m| + k_1 \right] |v_p - v_m|^2 \, d\Gamma. \quad (3.24) $$

Since $E$ is a cubic functional, local minima and local maxima can both occur depending on the signs and magnitudes of the constants $k_1$ and $k_2$. A local minimum $u \in V_{M_0}^{pm}$ is characterized by

$$ \forall v \in V_0^{pm} \quad dE(u; v) = 0 \quad \text{and} \quad \forall 0 \neq v \in V_0^{pm} \quad d^2 E(u; v; v) > 0 $$

and a local maximum $u \in V_{M_0}^{pm}$ by

$$ \forall v \in V_0^{pm} \quad dE(u; v) = 0 \quad \text{and} \quad \forall 0 \neq v \in V_0^{pm} \quad d^2 E(u; v; v) < 0. $$
Going back to the evolution equation (3.8) using the above conditions at the interface, we get the following system of equations

$$
\begin{align}
\frac{\partial c_p}{\partial t} &= \operatorname{div}(D_p \nabla c_p) \quad \text{in } \Omega_p, \quad \frac{\partial c_m}{\partial t} = \operatorname{div}(D_m \nabla c_m) \quad \text{in } \Omega_m \nonumber \\
c_p(x, 0) &= M_0 / |\Omega_p| \chi_{\Omega_p}(x) \quad \text{in } \Omega_p, \quad c_m(x, 0) = 0 \quad \text{in } \Omega_m \nonumber \\
D_p \frac{\partial c_p}{\partial n_p} + k_2 \frac{|\Omega_p|}{M_0} |c_p - c_m| (c_p - c_m) + k_1 (c_p - c_m) &= 0 \quad \text{on } \Gamma_{\text{int}} \nonumber \\
D_m \frac{\partial c_m}{\partial n_m} - \left[ k_2 \frac{|\Omega_p|}{M_0} |c_p - c_m| (c_p - c_m) + k_1 (c_p - c_m) \right] &= 0 \quad \text{on } \Gamma_{\text{int}} \tag{3.25} \\
D_p \frac{\partial c_p}{\partial n_p} &= 0 \quad \text{on } \Gamma_p \setminus \Gamma_{\text{int}}, \quad D_m \frac{\partial c_m}{\partial n_m} = 0 \quad \text{on } \Gamma_m \setminus \Gamma_{\text{int}} \nonumber \\
\int_{\Omega_p} c_p \, dx + \int_{\Omega_m} c_m \, dx &= M_0. \nonumber
\end{align}
$$

The nonlinear condition on $\Gamma_{\text{int}}$

$$
\begin{align*}
D_m \frac{\partial c_m}{\partial n_m} &= k_2 \frac{|\Omega_p|}{M_0} |c_p - c_m| (c_p - c_m) + k_1 (c_p - c_m) \\
&= \underbrace{\left( k_2 \frac{|\Omega_p|}{M_0} |c_p - c_m| + k_1 \right)}_{k(c)} (c_p - c_m)
\end{align*}
$$

says that $k(c)$ is an affine function of the size of the jump. This means that the rate of transfer of the product across the interface is large when the absolute value of the concentration jump is large. Assuming that $k_1 \ge 0$, when $k_2 > 0$ it decreases to $k_1$ when the size of the jump goes to zero; when $k_2 < 0$ it increases to $k_1$ when the size of the jump goes to zero.

**Remark 1.** When $k_2 > 0$, it would not be appropriate to remove the absolute value on $c_p - c_m$ in the term $k'$ of the previous identity. This would give the expression

$$
D_m \frac{\partial c_m}{\partial n_m} = \underbrace{\left( k_2 \frac{|\Omega_p|}{M_0} (c_p - c_m) + k_1 \right)}_{k'(c)} (c_p - c_m),
$$

where, if the size of the jump is large, $k'(c) > 0$ is large, $\partial c_m / \partial n_m > 0$ is large, and the diffusion of product would be from the medium to the polymer even when $c_p > c_m$, that is, when the concentration in the polymer is larger than the one in the medium. However, it is interesting to note that various behaviours can be modelled by replacing $|c_p - c_m|$ by the plus $[c_p - c_m]^+ = \max\{0, c_p - c_m\}$ or the minus $[c_p - c_m]^- = \max\{0, -(c_p - c_m)\}$ functions or introducing a threshold $\theta > 0 \max\{|c_p - c_m| - \tau, \theta\}$.
### 3.3 Relation between the PDE and the ODE models

Since $|\Omega_p|$ is much smaller than $|\Omega_m|$, this last equation is related to the quadratic ODE model by making the same assumptions on the concentrations on $\Gamma_{int}$ as in the previous section:

$$c_p(x, t) \approx \frac{1}{|\Omega_p|} \int_{\Omega_p} c_p(x, t) dx \quad \text{and} \quad c_m(x, t) \approx \frac{1}{|\Omega_m|} \int_{\Omega_m} c_m(x, t) dx \quad (3.26) \\ \Rightarrow c_p(x, t) - c_m(x, t) \approx \frac{1}{|\Omega_p|} [M_0 - M_m(t)],$$

where

$$M_m(t) \stackrel{\text{def}}{=} \int_{\Omega_m} c_m(x, t) dx \quad (3.27)$$

is the mass released at time $t$ in the medium and

$$\frac{dM_m}{dt}(t) = \frac{|\Gamma_{int}|}{|\Omega_p|} \left[ k_1 + \frac{k_2}{M_0} |M_0 - M_m(t)| \right] (M_0 - M_m(t)) \quad (3.28)$$

$$\Rightarrow \frac{dm_m}{dt}(t) = \frac{1}{h} [k_1 + k_2 |1 - m_m(t)|] (1 - m_m(t)), \quad m_m(t) \stackrel{\text{def}}{=} \frac{M_m(t)}{M_0}, \quad (3.29)$$

where $h = |\Omega_p|/|\Gamma_{int}|$ is the thickness of the polymer. This would correspond to $A_1 = k_1/h$ and $A_2 = k_2/h$ in the ODE model. The thickness $h$ is an important parameter: the thinner the polymer the faster the release. If $k_1$ and $k_2$ are constants, $m_m$ can be normalized through the change of variable $t \mapsto \tau = t/h$.

## References

1. R. P. Batycky, J. Hanes, R. Langer, D. A. Edwards, A theoretical model of erosion and macromolecular drug release from biodegrading microspheres, J. Pharm. Sci. 86 (1997) 1464–1477.

2. G. Blanchet, M. C. Delfour, and A. Garon, Quadratic models to fit experimental data of Paclitaxel release kinetics from biodegradable polymers, SIAM J. on Applied Mathematics (Special Issue on Mathematical Modeling of Controlled Drug Delivery) 71 (6) (2011), 2269–2286.

3. É. Bourgeois and M. C. Delfour, General patterns and asymptotic dose in the design of coated stents, Computer Methods in Biomechanics and Biomedical Engineering 11, (4) (2008), 323–334.

4. D. Cioranescu, A. Damlamian, G. Griso, D. Onofrei, The periodic unfolding method for perforated domains and Neumann sieve models, J. Math. Pures Appl. 89 (2008), 248–277

5. A. Damlamian, Le problème de la passoire de Neumann, (French) [The Neumann sieve problem], Rend. Sem. Mat. Univ. Politec. Torino 43 (1985), 427–450.

6. M. C. Delfour, Drug release kinetics from biodegradable polymers via partial differential equations models, Acta Appl Math 118 (2012), 161–183.

7. M. C. Delfour and A. Garon, New equations for the dose under pulsative/periodic conditions in the design of coated stents, Computer Methods in Biomechanics and Biomedical Engineering 13 (2010), (1), 19–34.
8. M. C. Delfour, A. Garon, and V. Longo, *Modeling and design of stents to optimize the effect of the dose*, SIAM J. on Applied Mathematics **65**, (3) (2005), 858–881.

9. M. C. Delfour, G. Payre, and J.-P. Zolésio, *Approximation of nonlinear problems associated with radiating bodies in space*, SIAM J. on Numerical Analysis **24** (1987), 1077–1094.

10. N. Faisant, J. Akiki, J. Siepmann, J. P. Benoit, J. Siepmann, *Effects of the type of release medium on drug release from PLGA-based microparticles: experiment and theory*, Int. J. Pharm. **314** (2006) 189–197.

11. A. Farb, P. F. Heller, S. Shroff, L. Cheng, F. D. Kolodgie, A. J. Carter, D. S. Scott, J. Froehlich, and R. Virmani, *Pathological analysis of local delivery of paclitaxel via a polymer-coated stent*, Circulation **104**, (4), (2001), 473–479.

12. A. Gopferich, *Polymer bulk erosion*, Macromolecules **30** (1997) 2598–2604.

13. T. Higuchi, *Mechanism of sustained action mediation: theoretical analysis of rate of release of solid drugs dispersed in solid matrices*, J. Pharm. Sci. **52** (1963) 1145–1149.

14. M. Joner, A. V. Finn, A. Farb, E. K. Mont, F. D. Kolodgie, E. Ladich, et al., *Pathology of drug-eluting stents in humans – delayed healing and late thrombotic risk*, J. Am. Coll. Cardiol. **48** (1) (2006), 193–202.

15. T. Kataoka, E. Grube, Y. Honda, Y. Morino, S.-H. Hur, H. N. Bonneau, A. Colombo, C. Di Mario, G. Guagliumi, K. E. Hauptmann, M. R. Pitney, A. J. Lansky, S. H. Stertzer, P. G. Yock, and P. J. Fitzgerald, *γ-Hexanoyltaxol-Eluting Stent for Prevention of Neointimal Growth: An Intravascular Ultrasound Analysis From the Study to Compare REstenosis rate between QueST and QuaDS-QP2 (SCORE)*, Circulation **106** (2002), 1788–1793.

16. L. L. Lao and S. S. Venkatraman, *Adjustable paclitaxel release kinetics and its efficacy to inhibit smooth muscle cells proliferation*, J. Control. Release **130** (2008), 9–14.

17. L. L. Lao, S. S. Venkatraman, *Paclitaxel release from single and double layered poly (DL-lactide-co-glycolide)/poly (L-lactide) film for biodegradable coronary stent application*, J. Biomed. Mater. Res. A Volume **87A** (1) (2008), 1–7.

18. L. L. Lao, S. S. Venkatraman, N. A. Peppas, *Modeling of drug release from biodegradable polymer blends*, Eur. J. Pharm. Biopharm. **70** (2008), 796–803.

19. L. L. Lao, S. S. Venkatraman, N. A. Peppas, *A novel model and experimental analysis of hydrophilic and hydrophobic agent release from biodegradable polymers*, J. Biomed. Mater. Res. A **90** (4) (2009), 1054–1065.

20. V. Lemaire, J. Bélair, P. Hildgen, *Structural modeling of a drug release from biodegradable porous matrices based on a combined diffusion/erosion process*, Int. J. Pharm. **258** (2003) 95–107.

21. E. Regar, G. Sianos, and P. W. Serruys, *Stent development and local drug delivery*, British Medical Bulletin **59** (1) (2001), 227–248.

22. J. Siepmann and A. Gopferich, *Mathematical modeling of bioerodable, polymeric drug delivery systems*, Adv. Drug Deliv. Rev. **48** (2001), 229–247.