# An Inverse Problem Solution Scheme for Solving the Optimization Problem of Drug-Controlled Release from Multilaminated Devices

Authors: Xinming Zhang
Journal: Computational and Mathematical Methods in Medicine
 (2020)

---

## Abstract

The optimization problem of drug release based on the multilaminated drug-controlled release devices has been solved in this paper under the inverse problem solution scheme. From the viewpoint of inverse problem, the solution of optimization problem can be regarded as the solution problem of a Fredholm integral equation of first kind. The solution of the Fredholm integral equation of first kind is a well-known ill-posed problem. In order to solve the severe ill-posedness, a modified regularization method is presented based on the Tikhonov regularization method and the truncated singular value decomposition method. The convergence analysis of the modified regularization method is also given. The optimization results of the initial drug concentration distribution obtained by the modified regularization method demonstrate that the inverse problem solution scheme proposed in this paper has the advantages of the numerical accuracy and antinoise property.


---

## Body

## 1. Introduction

It is well known that the controlled release device is usually used to regulate the release of active material for maintaining a preset concentration of the active material for a specified period of time. In recent years, the technique has been widely used in many fields, including drug, food, and cosmetics [1, 2], due to its safety and effectiveness. Especially, in the drug research field, the controlled release device has been of considerable interest to many scholars. Generally, the burst effect is not a desired delivery profile for drug release behavior, which will cause overdose of drug and systemic toxicity [3]. During the last three decades, several matrix geometries have been used to avoid this undesired effect, in which, the multilayer matrix device is an effective and simple device geometry[4–6]. For the multilayer matrix devices, the matrix core, containing the drug, is covered by one or more modulating layers that act as rate-controlling barriers. Researches have revealed that for the multilayer matrix devices, the reasonable initial drug loading can efficiently control the burst effect [7].For understanding and simulating the real release mechanisms from these complicated multilaminated matrix devices better, application of mathematical models is a good choice [8–10]. In 1961, Higuchi [11], who is regarded as the “father of mathematical modeling of drug release”, proposed a simple but surprising equation to quantify the drug release from monolithic dispersions with slab geometry: the famous Higuchi equation. After that, more and more scientists paid their attentions to the mathematical modeling of drug release, which can significantly contribute to the product development and mechanism understanding. Especially, over the last thirty years, there are many new progresses for different controlled release devices. Helbling et al. [12] used the refined integral method to study the controlled dispersed drug release from planar nonerodible polymeric matrices and derived an analytical solution in 2010 and then presented a novel mathematical model for drug-controlled release based on the torus-shaped single-layer devices in 2011. Streubel et al. [13] developed a new multilayer matrix tablets to obtain the bimodal drug release profiles and used theophylline and acetaminophen as the model drugs to investigate the effects of some parameters on the resulting release rates. Yin and Li [14] introduced the fractional calculus to model the anomalous diffusions and presented some new mathematical models to describe the drug release based on degradable and nondegradable slab matrices. Wu and Zhou [15] studied several factors affecting the kinetics of diffusional drug release based on the finite element method. Nevertheless, most of the previous methods focus on how to forecast the drug release profiles based on the given parameters. There have been few studies to investigate how to choose the suitable control parameters, in multilayer devices, to make the drug release profile as close to the desired release profile as possible. As far as we known, in 1998 and 1999, articles [5, 16] made relatively big contribution to this area firstly. In these two manuscripts, Lu and his coworkers employed a formal optimization approach to correctly determine the initial drug concentration in the layer so as to coincide with the required release profile as much as possible. To obtain the desired release rates in multilaminated drug delivery devices, Georgiadis and Kostoglou [7] presented a systematic optimization framework based on a simple mathematical model. Nauman et al. [17] designed a multilaminated drug delivery device, which has two or three layers with different initial drug concentration distributions. By adjusting the parameters in this device, all kinds of drug release profiles can be obtained. As shown in previous articles, the optimizations of the available control parameters have been performed in the frame of optimization. However, we can also treat this problem from another viewpoint, that is, the viewpoint of inverse problem. In fact, many practical problems in different fields can be reduced into the solution of inversion problems, such as geophysical exploration and medical imaging. In this article, from the viewpoint of the inverse problem, we adopt a solution scheme of inverse problem to solve the optimization problem of multilaminated drug-controlled release, which has been transformed into a diffusion equation initial value inverse problem. A classical regularization method, that is Tikhonov regularization method [18], and its variant have been attempted to solve the inverse problem. In 1963, Tikhonov proposed the Tikhonov regularization method firstly. Since then, many scholars turned their attentions to this rather effective method, developed many different variants in different space settings, and applied these regularization methods to solve all kinds of problems in science and engineering fields [19–22]. Fu et al. [23] applied the Fourier method to solve some ill-posed problems and systematically considered a posteriori choice of the regularization parameter; Zheng and Wei [24] used a spectral regularization method to solve the Cauchy problem of TFADE based on the solution given by the Fourier method; Zhao [25] introduced a mollification method to solve the ill-posed problem by using the expanded Hermit functions; Cheng et al. [26] presented an optimal filtering method to approximate a Cauchy problem for the Helmholtz equation in a rectangle and showed the Hölder-type error estimate; Bonesky et al. [27] applied an adaptive wavelet algorithms to solve an inverse parabolic problem describing the industrial process of melting iron ore in a steel furnace; Zhang and Li [28] established a new regularization method to solve the ill-posed problem based on the singular system theory of compact operator.In this article, we propose a modified regularization method based on a new regularizing filter function obtained by combining the method presented in the paper [28] with the truncated singular value decomposition method. We also show the convergence analysis of the proposed method. The regularization parameters are determined with the L-curve method suggested by Hansen in [29]. Then, we apply the modified regularization method to the optimization problem of the initial drug concentration distribution based on the multilaminated drug-controlled release devices and obtain some better results.In the reminder of the paper, it is organized as follows. We describe the mathematical model for the multilaminated drug-controlled release devices and the corresponding inverse problem solution scheme for the optimization of drug-controlled release in Section2. In Section 3, a modified regularization method is described in detail. This is followed by numerical simulation in Section 4 and Section 5, and lastly, the conclusion is indicated in Section 6.


## 2. Mathematical Model and Inverse Problem Solution Scheme

### 2.1. Mathematical Model

Figure1 shows a multilaminated drug-controlled release device, which has N layers. The device has a thickness L and initial drug concentration VX. It is sealed at the leftmost side by an impermeable layer and contacts with outside through the rightmost side. It is assumed that the device is not significantly swelling and eroding during drug release. Here, Cii=1,2,⋯,N are the drug concentrations of each layer, respectively. The present analysis focuses only on the case of low drug concentration.Mathematically, this problem can be modeled as one-dimensional partial differential equation based on Fick’s second law of diffusion:(1)∂C∂τ=∂∂XD∂C∂X,where C is the drug concentration, D is the diffusivity, and X and τ are the position and release time for one-dimensional diffusional processes, respectively.Figure 1
Drug release from a multilaminated drug-controlled release device.Under the assumption that zero flux and zero concentration are prescribed at the interface with impermeable layer and the environment, respectively, thus, the boundary conditions are as follows:(2)∂C∂XX=0=0,τ>0,(3)Cτ,L=0,τ>0.The initial conditions are imposed as follows:(4)C0,X=VX,τ=0,0<X<L.The flux of drug is also defined as follows:(5)Jτ,L=−DX=L∂C∂XX=L,τ>0,where VX denotes the initial drug concentration.


### 2.2. The Inverse Problem Solution Scheme

For the mathematical model (Equations (1)–(5)), if the initial conditions and the boundary conditions are known, the process to compute the concentration distribute functions Cτ,X is a forward problem, which is a well-posed problem. Conversely, how to identify the initial Condition (4) is a classical inverse problem, if we know the boundary Conditions (2) and (3) and additional Condition (5). Generally speaking, the inverse problem is ill-posed and always needs to be solved with some special algorithms, e.g., the Tikhonov regularization method.Assume that the diffusivity is constant, we use the following dimensionless processing to simplify the computation process:c=C/C0, v=VX/C0, x=X/L, t=D0τ/L2, j=JL/D0C0, d=D/D0=1, where C0 is a reference concentration and D0 is a reference diffusivity.Thus, the previous mathematical model (Equations (1)–(5)) can be rewritten as follows:
(6)∂c∂t=∂2c∂x2.Boundary conditions:(7)∂c∂xx=0=0,t>0,ct,1=0,t>0.Initial conditions:(8)c0,x=vx,t=0,0<x<1,Additional conditions:(9)jt,1=−∂c∂xx=1=jt,t>0.In fact, from the viewpoint of inverse problem, the problem to determine the initial conditionsvt based on the above mathematical model (Equations (6)–(9)) is a diffusion equation inverse problem and can further come down to a solution problem of Fredholm integral equation of first kind.The first step is to solve Equations (6)–(8), which is a diffusion equation initial boundary value problem. The method of separating variables leads to the analytical solution:
(10)ct,x=∑k=0∞2e−k+1/2π2tcosk+12πx∫01vxcosk+12πxdx.From Equation (10), the flux jt, namely, the additional conditions, in Equation (9), can be determined by differentiation with respect to x as follows:
(11)jt=−∂c∂xx=1=∑k=0∞2e−k+1/2π2tk+12πsink+12πx∫01vxcosk+12πxdx=∑k=0∞2−1kk+12πe−k+1/2π2t∫01vxcosk+12πxdx=2∫01∑k=0∞−1kk+12πe−k+1/2π2tcosk+12πxvxdx.Then, using Equation (11), the previous inverse problem (Equations (6)–(9)) for determining initial conditions can further be transformed into the following Fredholm integral equation of first kind:
(12)2∫01∑k=0∞−1kk+12πe−k+1/2π2tcosk+12πxvxdx=jt,where vx is the unknown function and the kernel function is
(13)kx,t=∑k=0∞−1kk+12πe−k+1/2π2tcosk+12πx.Thus, successful solution of the Fredholm integral equation of first kind will lead to the effective identification of the initial conditions in the mathematical model (Equations (6)–(9)), that is, the optimization of the initial drug concentration in the optimization problem of drug release based on the multilaminated drug-controlled release devices. However, the solution of the Fredholm integral equation of first kind is a well-known ill-posed problem, especially, for the case of input data with noise. We usually need a suitable and specific method to solve it. In the next section, we propose a new regularization method to solve Equation (12).


## 2.1. Mathematical Model

Figure1 shows a multilaminated drug-controlled release device, which has N layers. The device has a thickness L and initial drug concentration VX. It is sealed at the leftmost side by an impermeable layer and contacts with outside through the rightmost side. It is assumed that the device is not significantly swelling and eroding during drug release. Here, Cii=1,2,⋯,N are the drug concentrations of each layer, respectively. The present analysis focuses only on the case of low drug concentration.Mathematically, this problem can be modeled as one-dimensional partial differential equation based on Fick’s second law of diffusion:(1)∂C∂τ=∂∂XD∂C∂X,where C is the drug concentration, D is the diffusivity, and X and τ are the position and release time for one-dimensional diffusional processes, respectively.Figure 1
Drug release from a multilaminated drug-controlled release device.Under the assumption that zero flux and zero concentration are prescribed at the interface with impermeable layer and the environment, respectively, thus, the boundary conditions are as follows:(2)∂C∂XX=0=0,τ>0,(3)Cτ,L=0,τ>0.The initial conditions are imposed as follows:(4)C0,X=VX,τ=0,0<X<L.The flux of drug is also defined as follows:(5)Jτ,L=−DX=L∂C∂XX=L,τ>0,where VX denotes the initial drug concentration.


## 2.2. The Inverse Problem Solution Scheme

For the mathematical model (Equations (1)–(5)), if the initial conditions and the boundary conditions are known, the process to compute the concentration distribute functions Cτ,X is a forward problem, which is a well-posed problem. Conversely, how to identify the initial Condition (4) is a classical inverse problem, if we know the boundary Conditions (2) and (3) and additional Condition (5). Generally speaking, the inverse problem is ill-posed and always needs to be solved with some special algorithms, e.g., the Tikhonov regularization method.Assume that the diffusivity is constant, we use the following dimensionless processing to simplify the computation process:c=C/C0, v=VX/C0, x=X/L, t=D0τ/L2, j=JL/D0C0, d=D/D0=1, where C0 is a reference concentration and D0 is a reference diffusivity.Thus, the previous mathematical model (Equations (1)–(5)) can be rewritten as follows:
(6)∂c∂t=∂2c∂x2.Boundary conditions:(7)∂c∂xx=0=0,t>0,ct,1=0,t>0.Initial conditions:(8)c0,x=vx,t=0,0<x<1,Additional conditions:(9)jt,1=−∂c∂xx=1=jt,t>0.In fact, from the viewpoint of inverse problem, the problem to determine the initial conditionsvt based on the above mathematical model (Equations (6)–(9)) is a diffusion equation inverse problem and can further come down to a solution problem of Fredholm integral equation of first kind.The first step is to solve Equations (6)–(8), which is a diffusion equation initial boundary value problem. The method of separating variables leads to the analytical solution:
(10)ct,x=∑k=0∞2e−k+1/2π2tcosk+12πx∫01vxcosk+12πxdx.From Equation (10), the flux jt, namely, the additional conditions, in Equation (9), can be determined by differentiation with respect to x as follows:
(11)jt=−∂c∂xx=1=∑k=0∞2e−k+1/2π2tk+12πsink+12πx∫01vxcosk+12πxdx=∑k=0∞2−1kk+12πe−k+1/2π2t∫01vxcosk+12πxdx=2∫01∑k=0∞−1kk+12πe−k+1/2π2tcosk+12πxvxdx.Then, using Equation (11), the previous inverse problem (Equations (6)–(9)) for determining initial conditions can further be transformed into the following Fredholm integral equation of first kind:
(12)2∫01∑k=0∞−1kk+12πe−k+1/2π2tcosk+12πxvxdx=jt,where vx is the unknown function and the kernel function is
(13)kx,t=∑k=0∞−1kk+12πe−k+1/2π2tcosk+12πx.Thus, successful solution of the Fredholm integral equation of first kind will lead to the effective identification of the initial conditions in the mathematical model (Equations (6)–(9)), that is, the optimization of the initial drug concentration in the optimization problem of drug release based on the multilaminated drug-controlled release devices. However, the solution of the Fredholm integral equation of first kind is a well-known ill-posed problem, especially, for the case of input data with noise. We usually need a suitable and specific method to solve it. In the next section, we propose a new regularization method to solve Equation (12).


## 3. The Modified Regularization Method

Kirsch [30] proposed the concept of regularizing filter function to investigate the ill-posed problem. As shown in [30], for the regularizing filter function, we have the following Lemma:Lemma 1.
LetX and Y be both Hilbert spaces, K:X→Y be compact with singular system μi,xi,yi, and q:0,+∞×0,K→R be a function with the following properties:
(1)
∣qα,μ∣≤1, ∀α∈0,+∞,∀μ∈0,K(2)
For∀α∈0,+∞, there exists cα>0, such that ∣qα,μ∣≤cαμ,∀μ∈0,K(3)
limα→0qα,μ=1, ∀μ∈0,K
Then, the operatorRα:Y→Xis defined by
(14)Rαy=∑i=1∞qα,μiμiy,yixi,which is a regularization strategy, where Rα≤cα. We call the function qα,μ satisfying the previous three properties regularizing filter function for K.Based on the singular system theory of the compact operator, the regularized solution of the famous Tikhonov regularization method can be expressed in the following formula:(15)xαδ=Rαyδ=∑i=1∞qα,μiμiyδ,yixi,where qα,μ=μ2/α+μ2, (α>0,0<μ≤K) is the Tikhonov regularizing filter function, and yδ is the disturbed data.In the paper [28], Zhang and Li presented an improved regularization method based on a new regularizing filter function qα,μ=μσ/α+μσr1/r,r>0,∈σ≥1 and proved that the regularization solution can achieve the optimal asymptotic convergence rate by selecting reasonable regularization parameters.Similarly, in this paper, a new regularization method is introduced based on the previous improved regularization method and the truncated singular value decomposition (TSVD) regularization method. The convergence and the optimal asymptotic order of the new regularized solution are also obtained.


### 3.1. A New Regularizing Filter Function

The regularizing filter function corresponding to the improved regularization method proposed in literature [28] is as follows:
(16)qα,μ=μσα+μσr1/r,r>0,σ≥1.The filter function of the truncated singular value decomposition regularization method is as follows:(17)qα,μ=1μ≥α0μ<α.Based on the previous two methods, we introduce a new regularizing filter function, which is defined as follows:(18)qα,μ=1μσr≥αμσα+μσr1/rμσr<α,where α>0,0<μ≤K,r>0,σ≥1.Using Equation (18), we can make the large singular value not be modified and the small singular value be filtered. Thus, the modified regularization method adopting the new regularizing filter function will lead to a more accurate regularized solution.For the new regularizing filter function, we prove the following theorem.Theorem 1.
The functionqα,μ defined by Equation (18) satisfies the assumptions (1), (2), and (3) of Lemma1, namely, the function qα,μ defined by Equation (18) is a regularizing filter function.Proof.
(1)
It is sufficient to consider the caseμσr<α. In this case, qα,μ=μσ/α+μσr1/r<1 because μσ=μσr1/r<α+μσr1/r(2)
For the case ofμσr≥α, we conclude that μσr/α≥1, that is, μ/α1/σr≥1; this proves that qα,μ=1≤1/α1/σrμ
For the case ofμσr<α , we have qα,μ=μσ/α+μσr1/r.
The choiceq=σ/σ−1 leads to
(19)qα,μ=μα+μr1/r≤μα1/r.
Ifσ>1, let p=σ and q=σ/σ−1. Thanks to the Yong inequality α+μσr≥α1/p⋅μσr/q, we have α+μσr≥α1/σ⋅μrσ−1, which yields
(20)μσα+μσr1/r≤μσα1/σr⋅μσ−1=μα1/σr.
Thus, the inequalityqα,μ≤μ/α1/σr holds for ∀σ≥1. That is the property (2) in Lemma 1 holds withcα=1/α1/σr, for∀α>0.
(3)
It is obvious thatqα,μ=1, as α→0So, fromLemma1, we can know that the function qα,μ defined by Equation (18) is a regularizing filter function and the corresponding regularization operator Rα:Y→X is as follows:
(21)Rαy=∑i=1∞qα,μiμiy,yixi.Then, a modified regularization method can be constructed based on the regularization operatorRα from Equation (21), which can also be used to solve the ill-posed problem for the case of input data with noise. The regularized solution xαδ is therefore defined by xαδ=Rαyδ=∑i=1∞qα,μi/μiyδ,yixi.We combine the singular system theory of the compact operator with the previous theorem and show the following result for the regularized solution.


### 3.2. Error Analysis of the Regularized Solution

Theorem 2.
Letx+=K∗Kvz∈RK∗Kv, z∈X, with z≤E; for the regularization operators Rα:Y→X from Equation (21), we choose the regularized parameters αδ=cδ/Eσr/2v+1 for some c>0; then, the following estimate holds:
(22)xαδδ−x+=Oδ2v/2v+1,where x+ is the solution of Kx=y, xαδδ≔Rαyδ is the approximation of x, and K∗ denotes the adjoint operator of K.Proof.
The error betweenx and xαδδ is
(23)xαδ−x+≤Rα⋅δ+Rαy−x+,Theorem1 yields Rα≤cα=1/α1/σr. From RαKx=∑i=1∞qα,μi/μiKx,yixi, x=∑i=1∞x,xixi, and Kx,yi=x,K∗yi=μix,xi, we conclude that
(24)Rαy−x+2=∑i=1∞qα,μi−12⋅x+,xi2=∑i=1∞qα,μi−12⋅K∗Kvz,xi2=∑i=1∞qα,μi−12⋅∑j=1∞μj2vz,xjxj,xi2=∑i=1∞qα,μi−12⋅μi2vz,xi2=∑i=1∞qα,μi−12⋅μi4v⋅z,xi2.Ifμσr<α, qα,μi−1<1 because 0<qα,μ<1; thus,
(25)qα,μi−1⋅μi2v<μi2v=μiσr2v/σr<α2v/σr,Ifμσr≥α, qα,μ=1, so
(26)qα,μi−1⋅μi2v=0⋅μi2v=0<α2v/σr.Then,qα,μi−1⋅μi2v<α2v/σr holds for all two cases.Therefore,(27)Rαy−x+2<α4v/σr∑i=1∞z,xi2=α4v/σr⋅z2≤α4v/σr⋅E2,that is Rαy−x+<α2v/σr⋅E. Thus, we have shown that
(28)xαδ−x+≤1α1/σr⋅δ+α2v/σr⋅E.The choice ofαδ=cδ/Eσr/2v+1 leads to the corresponding estimate
(29)xαδδ−x+≤δ⋅cδEσr/2v+1−1/σr+cδEσr/2v+12v/σr⋅E=δ⋅c−1/σr⋅δE−1/2v+1+c2v/σr⋅δE2v/2v+1⋅E=c−1/σr⋅E1/2v+1⋅δ1−1/2v+1+c2v/σr⋅δ2v/2v+1⋅E1−2v/2v+1=c−1/σr+c2v/σr⋅E1/2v+1⋅δ2v/2v+1.This yields(30)xαδδ−x+=Oδ2v/2v+1.


## 3.1. A New Regularizing Filter Function

The regularizing filter function corresponding to the improved regularization method proposed in literature [28] is as follows:
(16)qα,μ=μσα+μσr1/r,r>0,σ≥1.The filter function of the truncated singular value decomposition regularization method is as follows:(17)qα,μ=1μ≥α0μ<α.Based on the previous two methods, we introduce a new regularizing filter function, which is defined as follows:(18)qα,μ=1μσr≥αμσα+μσr1/rμσr<α,where α>0,0<μ≤K,r>0,σ≥1.Using Equation (18), we can make the large singular value not be modified and the small singular value be filtered. Thus, the modified regularization method adopting the new regularizing filter function will lead to a more accurate regularized solution.For the new regularizing filter function, we prove the following theorem.Theorem 1.
The functionqα,μ defined by Equation (18) satisfies the assumptions (1), (2), and (3) of Lemma1, namely, the function qα,μ defined by Equation (18) is a regularizing filter function.Proof.
(1)
It is sufficient to consider the caseμσr<α. In this case, qα,μ=μσ/α+μσr1/r<1 because μσ=μσr1/r<α+μσr1/r(2)
For the case ofμσr≥α, we conclude that μσr/α≥1, that is, μ/α1/σr≥1; this proves that qα,μ=1≤1/α1/σrμ
For the case ofμσr<α , we have qα,μ=μσ/α+μσr1/r.
The choiceq=σ/σ−1 leads to
(19)qα,μ=μα+μr1/r≤μα1/r.
Ifσ>1, let p=σ and q=σ/σ−1. Thanks to the Yong inequality α+μσr≥α1/p⋅μσr/q, we have α+μσr≥α1/σ⋅μrσ−1, which yields
(20)μσα+μσr1/r≤μσα1/σr⋅μσ−1=μα1/σr.
Thus, the inequalityqα,μ≤μ/α1/σr holds for ∀σ≥1. That is the property (2) in Lemma 1 holds withcα=1/α1/σr, for∀α>0.
(3)
It is obvious thatqα,μ=1, as α→0So, fromLemma1, we can know that the function qα,μ defined by Equation (18) is a regularizing filter function and the corresponding regularization operator Rα:Y→X is as follows:
(21)Rαy=∑i=1∞qα,μiμiy,yixi.Then, a modified regularization method can be constructed based on the regularization operatorRα from Equation (21), which can also be used to solve the ill-posed problem for the case of input data with noise. The regularized solution xαδ is therefore defined by xαδ=Rαyδ=∑i=1∞qα,μi/μiyδ,yixi.We combine the singular system theory of the compact operator with the previous theorem and show the following result for the regularized solution.


## 3.2. Error Analysis of the Regularized Solution

Theorem 2.
Letx+=K∗Kvz∈RK∗Kv, z∈X, with z≤E; for the regularization operators Rα:Y→X from Equation (21), we choose the regularized parameters αδ=cδ/Eσr/2v+1 for some c>0; then, the following estimate holds:
(22)xαδδ−x+=Oδ2v/2v+1,where x+ is the solution of Kx=y, xαδδ≔Rαyδ is the approximation of x, and K∗ denotes the adjoint operator of K.Proof.
The error betweenx and xαδδ is
(23)xαδ−x+≤Rα⋅δ+Rαy−x+,Theorem1 yields Rα≤cα=1/α1/σr. From RαKx=∑i=1∞qα,μi/μiKx,yixi, x=∑i=1∞x,xixi, and Kx,yi=x,K∗yi=μix,xi, we conclude that
(24)Rαy−x+2=∑i=1∞qα,μi−12⋅x+,xi2=∑i=1∞qα,μi−12⋅K∗Kvz,xi2=∑i=1∞qα,μi−12⋅∑j=1∞μj2vz,xjxj,xi2=∑i=1∞qα,μi−12⋅μi2vz,xi2=∑i=1∞qα,μi−12⋅μi4v⋅z,xi2.Ifμσr<α, qα,μi−1<1 because 0<qα,μ<1; thus,
(25)qα,μi−1⋅μi2v<μi2v=μiσr2v/σr<α2v/σr,Ifμσr≥α, qα,μ=1, so
(26)qα,μi−1⋅μi2v=0⋅μi2v=0<α2v/σr.Then,qα,μi−1⋅μi2v<α2v/σr holds for all two cases.Therefore,(27)Rαy−x+2<α4v/σr∑i=1∞z,xi2=α4v/σr⋅z2≤α4v/σr⋅E2,that is Rαy−x+<α2v/σr⋅E. Thus, we have shown that
(28)xαδ−x+≤1α1/σr⋅δ+α2v/σr⋅E.The choice ofαδ=cδ/Eσr/2v+1 leads to the corresponding estimate
(29)xαδδ−x+≤δ⋅cδEσr/2v+1−1/σr+cδEσr/2v+12v/σr⋅E=δ⋅c−1/σr⋅δE−1/2v+1+c2v/σr⋅δE2v/2v+1⋅E=c−1/σr⋅E1/2v+1⋅δ1−1/2v+1+c2v/σr⋅δ2v/2v+1⋅E1−2v/2v+1=c−1/σr+c2v/σr⋅E1/2v+1⋅δ2v/2v+1.This yields(30)xαδδ−x+=Oδ2v/2v+1.


## 4. Solution of the Fredholm Integral Equation of First Kind

For validation purpose, in this section, the modified regularization method proposed in this paper is applied to solve the Fredholm integral equation of first kind.Example 1.
We first consider the following integral equation of first kind.(31)∫011+tsetsxsds=yt,0≤t≤1,
with the right-hand side and the analysis solution given byyt=et and xt=1.We discretize the integral equation by the compound trapezoidal formula and obtain the linear system.(32)AX=b,where A∈R100×100 and the error-free right-hand side b=yt0,⋯,ytmT∈R100. The associated contaminated vector b~ is given by the following formula.
(33)b~=b1+δ×rand,where δ denotes the noise level and rand is a number randomly generated within the interval 0,1. In this example, the noise level is assumed to be 0.001.We solve equationAX=b~ with the Tikhonov regularization method and our method, respectively. The regularization parameters are chosen by using the L-curve method. The corresponding L-curve for the Tikhonov regularization method is shown in Figure 2. The corner of L-curve is located at the points (AX−b~2, X2) with the regularization parameter α=3.7866×10−3. We set the parameters of the modified regularization method σ=4 and r=1.5 as suggested in the paper [28]. Figure 3 depicts the comparison between the results obtained by the classical Tikhonov regularization method and the exact solution. Figure 4 shows the corresponding L-curve for the modified regularization method, and the regularization parameter is 2.1147∗e−3. And we compare the results obtained by the modified regularization method and the exact solution in Figure 5. As we can see that, from these two results, the modified regularization method is more effective than the classical Tikhonov regularization method.Example 2.
Consider the following Fredholm integral equation of first kind.(34)∫01etsxsds=yt,0≤t≤1,which has a unique exact solution xt=e2t, and the right-hand side is given by yt=et+2−1/t+2.Figure 2
Regularization parameter choice (L-curve).Figure 3
The results of Tikhonov regularization method.Figure 4
Regularization parameter choice (L-curve).Figure 5
The results of modified regularization method.The parameters involved in the numerical simulation are the same as the previous example. We also apply the Tikhonov regularization method and our method to solve Equation (34), respectively. The L- curve for Tikhonov regularization method is shown in Figure 6. The corner of L-curve is located at the point for α=4.2365×10−3. Figure 7 depicts the comparison between the results obtained by the classical Tikhonov regularization method and the exact solution. Figure 8 shows the corresponding L-curve for the modified regularization method, and the regularization parameter is 2.7659∗e−3. Figure 9 gives the comparison between the results obtained by the modified regularization method and the exact solution. From Figures 7 and 9, we can conclude that the modified regularization method works better for this problem.Figure 6
Regularization parameter choice (L-curve).Figure 7
The results of Tikhonov regularization method.Figure 8
Regularization parameter choice (L-curve).Figure 9
The results of modified regularization method.From the previous numerical results, we can conclude that the modified regularization method is effective for the Fredholm integral equation of first kind for the case of input data with noise. Meanwhile, as shown in Section2, the determination of the initial drug concentration in the optimization problem of drug release based on the multilaminated drug-controlled release devices can be transformed into the solution of the Fredholm integral equation of first kind. So, in the next section, we adopt the proposed method to deal with the optimization problem of drug-controlled release from multilaminated devices.


## 5. Optimization of Drug-Controlled Release from Multilaminated Devices

### 5.1. The Optimization of Initial Drug Concentration

The initial drug concentration is an essential parameter in the multilaminated controlled release system, which can affect the drug release greatly. A reasonable initial drug concentration can lead to the drug release with the desired flux. Based on the inverse problem solution scheme, the initial drug concentration can be inverted from the known drug release flux. In the following, we will determine the initial drug concentration for three different cases with the proposed inverse problem solution scheme based on the Tikhonov regularization method (TRM) and the modified regularization method (MRM), respectively. The three different desired release profiles are shown in Figure10.Figure 10
Three different desired flux. Case 1: constant release rate. Case 2: linearly decreasing release rate. Case 3: nonlinear release rate.


#### 5.1.1. Case 1

We first consider a typical case. In this case, suppose that the desired flux is constant (jt=1,0≤t≤0.5). For the Fredholm integral equation of the first kind (Equation (12)), set the right-hand side equal to 1, that is, jt=1. The Tikhonov regularization method and modified regularization method are applied to solve this ill-posed problem. The inverse results are shown in Figure 11, and the computational drug release flux based on the inversed initial drug concentration are depicted in Figure 12. It is seen that from Figure 12, the optimal release profile obtained with TRM remains flatter pattern at initial stage and has a smaller deviation from ideal case. However, bigger error appears as the time increases. The mean square deviation for TRM can reach 0.2174. We can also observe that, for MRM, although the optimized release fluctuates at the initial stage, the deviation from the desired release remains smaller over the entire computation time. Compared with TRM, the mean square deviation for MRM is also relatively less, which is 0.1692.Figure 11
The inverse results. (a) Results obtained by TRM. (b) Results obtained by MRM.
(a)(b)Figure 12
The computational drug release flux based on the inversed initial drug concentration.


#### 5.1.2. Case 2

In some cases, we maybe desire to obtain an approximately linearly increasing profile, e.g., to build up a tolerance for the chemical material transmitted. In the following, the ideal release profile is given with the functionjt=1.5−2t,0≤t≤0.5. The inversed results obtained by using the previous two different regularization methods are shown in Figure 13 and the optimized release profiles are depicted in Figure 14. We can conclude that MRM has obvious advantages over the TRM. In Figure 14, the optimized release profile almost coincides with the desired release flux for MRM, whereas the performance of TRM remains poor. The mean square deviations of two methods are 0.1285 and 0.0164, respectively.Figure 13
The inverse results. (a) Results obtained by TRM. (b) Results obtained by MRM.
(a)(b)Figure 14
The optimized drug release flux based on the inversed initial drug concentration.


#### 5.1.3. Case 3

Some situations demand a nonlinearly release rate, e.g., linearly increasing follow by a constant release, and without burst. A typical example is the delivery of the some anticancer drug.For this case, the release rate function is(35)jt=24t0≤t≤0.051.20.05<t≤0.5.The inversed results obtained by using two different regularization methods are shown in Figure15 and the optimized release profiles are depicted in Figure 16. From these two figures, we can see that the performance of MRM is better than that of TRM. The mean square deviations of two methods are 0.3019 and 0.2164, respectively.Figure 15
The inverse results. (a) Results obtained by TRM. (b) Results obtained by MRM.
(a)(b)Figure 16
The optimized drug release flux based on the inversed initial drug concentration.


### 5.2. Antinoise Property Analysis

In practice, we cannot guarantee that the right-hand side of the Fredholm integral equation of the first kind (Equation (12)) is known exactly. Instead, jt usually contains some error, say, δ>0. Therefore, we will consider the initial drug concentration optimization problem for the case of the right-hand side with a perturbed data.Assume that we knowδ>0 and with perturbed data jδt satisfying jδt−jt≤δ. It is our aim to solve the following perturbed equation.
(36)2∫01∑k=0∞−1kk+12πe−k+1/2π2tcosk+12πxvδxdx=jδt.We use the modified regularization method to solve the above perturbed equation for three different cases. Figure17 shows the error between the inverted results with noise and that without noise against the different noise levels. We can conclude from Figure 17 that the error increases apparently with the increment of noise level. However, even for the noise level of δ=0.3, the results are acceptable. To further show the optimization results intuitively, we list the initial drug concentration optimization results with δ=0.1 for three different cases in Figures 18–20. These three pictures demonstrate that the modified regularization method we proposed can still work well. It means that the MRM has a good antinoise property.Figure 17
The error against the noise level.Figure 18
The inverse initial drug concentration for Case 1 withδ=0.1.Figure 19
The inverse initial drug concentration for Case 2 withδ=0.1.Figure 20
The inverse initial drug concentration for Case 3 withδ=0.1.


## 5.1. The Optimization of Initial Drug Concentration

The initial drug concentration is an essential parameter in the multilaminated controlled release system, which can affect the drug release greatly. A reasonable initial drug concentration can lead to the drug release with the desired flux. Based on the inverse problem solution scheme, the initial drug concentration can be inverted from the known drug release flux. In the following, we will determine the initial drug concentration for three different cases with the proposed inverse problem solution scheme based on the Tikhonov regularization method (TRM) and the modified regularization method (MRM), respectively. The three different desired release profiles are shown in Figure10.Figure 10
Three different desired flux. Case 1: constant release rate. Case 2: linearly decreasing release rate. Case 3: nonlinear release rate.


### 5.1.1. Case 1

We first consider a typical case. In this case, suppose that the desired flux is constant (jt=1,0≤t≤0.5). For the Fredholm integral equation of the first kind (Equation (12)), set the right-hand side equal to 1, that is, jt=1. The Tikhonov regularization method and modified regularization method are applied to solve this ill-posed problem. The inverse results are shown in Figure 11, and the computational drug release flux based on the inversed initial drug concentration are depicted in Figure 12. It is seen that from Figure 12, the optimal release profile obtained with TRM remains flatter pattern at initial stage and has a smaller deviation from ideal case. However, bigger error appears as the time increases. The mean square deviation for TRM can reach 0.2174. We can also observe that, for MRM, although the optimized release fluctuates at the initial stage, the deviation from the desired release remains smaller over the entire computation time. Compared with TRM, the mean square deviation for MRM is also relatively less, which is 0.1692.Figure 11
The inverse results. (a) Results obtained by TRM. (b) Results obtained by MRM.
(a)(b)Figure 12
The computational drug release flux based on the inversed initial drug concentration.


### 5.1.2. Case 2

In some cases, we maybe desire to obtain an approximately linearly increasing profile, e.g., to build up a tolerance for the chemical material transmitted. In the following, the ideal release profile is given with the functionjt=1.5−2t,0≤t≤0.5. The inversed results obtained by using the previous two different regularization methods are shown in Figure 13 and the optimized release profiles are depicted in Figure 14. We can conclude that MRM has obvious advantages over the TRM. In Figure 14, the optimized release profile almost coincides with the desired release flux for MRM, whereas the performance of TRM remains poor. The mean square deviations of two methods are 0.1285 and 0.0164, respectively.Figure 13
The inverse results. (a) Results obtained by TRM. (b) Results obtained by MRM.
(a)(b)Figure 14
The optimized drug release flux based on the inversed initial drug concentration.


### 5.1.3. Case 3

Some situations demand a nonlinearly release rate, e.g., linearly increasing follow by a constant release, and without burst. A typical example is the delivery of the some anticancer drug.For this case, the release rate function is(35)jt=24t0≤t≤0.051.20.05<t≤0.5.The inversed results obtained by using two different regularization methods are shown in Figure15 and the optimized release profiles are depicted in Figure 16. From these two figures, we can see that the performance of MRM is better than that of TRM. The mean square deviations of two methods are 0.3019 and 0.2164, respectively.Figure 15
The inverse results. (a) Results obtained by TRM. (b) Results obtained by MRM.
(a)(b)Figure 16
The optimized drug release flux based on the inversed initial drug concentration.


## 5.1.1. Case 1

We first consider a typical case. In this case, suppose that the desired flux is constant (jt=1,0≤t≤0.5). For the Fredholm integral equation of the first kind (Equation (12)), set the right-hand side equal to 1, that is, jt=1. The Tikhonov regularization method and modified regularization method are applied to solve this ill-posed problem. The inverse results are shown in Figure 11, and the computational drug release flux based on the inversed initial drug concentration are depicted in Figure 12. It is seen that from Figure 12, the optimal release profile obtained with TRM remains flatter pattern at initial stage and has a smaller deviation from ideal case. However, bigger error appears as the time increases. The mean square deviation for TRM can reach 0.2174. We can also observe that, for MRM, although the optimized release fluctuates at the initial stage, the deviation from the desired release remains smaller over the entire computation time. Compared with TRM, the mean square deviation for MRM is also relatively less, which is 0.1692.Figure 11
The inverse results. (a) Results obtained by TRM. (b) Results obtained by MRM.
(a)(b)Figure 12
The computational drug release flux based on the inversed initial drug concentration.


## 5.1.2. Case 2

In some cases, we maybe desire to obtain an approximately linearly increasing profile, e.g., to build up a tolerance for the chemical material transmitted. In the following, the ideal release profile is given with the functionjt=1.5−2t,0≤t≤0.5. The inversed results obtained by using the previous two different regularization methods are shown in Figure 13 and the optimized release profiles are depicted in Figure 14. We can conclude that MRM has obvious advantages over the TRM. In Figure 14, the optimized release profile almost coincides with the desired release flux for MRM, whereas the performance of TRM remains poor. The mean square deviations of two methods are 0.1285 and 0.0164, respectively.Figure 13
The inverse results. (a) Results obtained by TRM. (b) Results obtained by MRM.
(a)(b)Figure 14
The optimized drug release flux based on the inversed initial drug concentration.


## 5.1.3. Case 3

Some situations demand a nonlinearly release rate, e.g., linearly increasing follow by a constant release, and without burst. A typical example is the delivery of the some anticancer drug.For this case, the release rate function is(35)jt=24t0≤t≤0.051.20.05<t≤0.5.The inversed results obtained by using two different regularization methods are shown in Figure15 and the optimized release profiles are depicted in Figure 16. From these two figures, we can see that the performance of MRM is better than that of TRM. The mean square deviations of two methods are 0.3019 and 0.2164, respectively.Figure 15
The inverse results. (a) Results obtained by TRM. (b) Results obtained by MRM.
(a)(b)Figure 16
The optimized drug release flux based on the inversed initial drug concentration.


## 5.2. Antinoise Property Analysis

In practice, we cannot guarantee that the right-hand side of the Fredholm integral equation of the first kind (Equation (12)) is known exactly. Instead, jt usually contains some error, say, δ>0. Therefore, we will consider the initial drug concentration optimization problem for the case of the right-hand side with a perturbed data.Assume that we knowδ>0 and with perturbed data jδt satisfying jδt−jt≤δ. It is our aim to solve the following perturbed equation.
(36)2∫01∑k=0∞−1kk+12πe−k+1/2π2tcosk+12πxvδxdx=jδt.We use the modified regularization method to solve the above perturbed equation for three different cases. Figure17 shows the error between the inverted results with noise and that without noise against the different noise levels. We can conclude from Figure 17 that the error increases apparently with the increment of noise level. However, even for the noise level of δ=0.3, the results are acceptable. To further show the optimization results intuitively, we list the initial drug concentration optimization results with δ=0.1 for three different cases in Figures 18–20. These three pictures demonstrate that the modified regularization method we proposed can still work well. It means that the MRM has a good antinoise property.Figure 17
The error against the noise level.Figure 18
The inverse initial drug concentration for Case 1 withδ=0.1.Figure 19
The inverse initial drug concentration for Case 2 withδ=0.1.Figure 20
The inverse initial drug concentration for Case 3 withδ=0.1.


## 6. Conclusion

We have proposed a new viewpoint to solve the optimization problem of drug release based on the multilaminated drug-controlled release device, that is, inverse problem solution scheme. The objective of this paper is to show that the inverse problem solution scheme is effective for the optimization problem of drug release. Based on the inverse problem solution scheme, the optimization problem of drug release can be transformed into the diffusion equation initial value inverse problem and further converted to the Fredholm integral equation of first kind. The solution of Fredholm integral equation of first kind is an ill-posed problem, which have to be solved by suitable regularization method.To solve this ill-posed problem, we introduce a new regularizing filter function and propose a modified regularization method. The error analysis of the regularized solution obtained by the proposed method is also verified. Furthermore, for three various desired release flux, the modified regularization method is applied to inverse the initial drug concentration. As seen in the examples, the method proposed in this paper has been successful at inverting the initial drug concentration. This demonstrates that the modified regularization approach is well suited to solving this ill-posed problem.Also shown in this paper is the result that the modified regularization method has a better antinoise property for the initial drug concentration estimation. With 10% noise, the results obtained with the MRM are satisfactory.In all, the better results obtained in this paper mean that the inverse problem solution scheme exhibits its effectiveness and superiority, for the optimization problem of drug release based on the multilaminated drug-controlled release device, to some extent in both theoretical research and numerical simulation. There is a good potential that the proposed method can be employed to solve more complicated cases, such as multiparameter identification and high-dimensional problem. And this is an important direction for us to face in future.



---

*Source: 8380691-2020-08-01.xml*