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

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.


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][5][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][9][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][20][21][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 Section 2. 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.

Mathematical Model and Inverse Problem
Solution Scheme 2.1. Mathematical Model. Figure 1 shows a multilaminated drug-controlled release device, which has N layers. The device has a thickness L and initial drug concentration VðXÞ. 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, C i ði = 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 onedimensional partial differential equation based on Fick's second law of diffusion: 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. 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: The initial conditions are imposed as follows: The flux of drug is also defined as follows: where VðXÞ denotes the initial drug concentration.  (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: Thus, the previous mathematical model (Equations (1)-(5)) can be rewritten as follows: Boundary conditions: Initial conditions: Additional conditions: In fact, from the viewpoint of inverse problem, the problem to determine the initial conditions vðtÞ 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: From Equation (10), the flux jðtÞ, namely, the additional conditions, in Equation (9), can be determined by differentiation with respect to x as follows: 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: where vðxÞ is the unknown function and the kernel function is 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).

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: Figure 1: Drug release from a multilaminated drug-controlled release device.
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: where qðα, μÞ = μ 2 /ðα + μ 2 Þ, (α > 0, 0 < μ ≤ kKk) 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ðα, μÞ = μ σ /ðα + μ σr Þ 1/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.

A New Regularizing Filter Function.
The regularizing filter function corresponding to the improved regularization method proposed in literature [28] is as follows: The filter function of the truncated singular value decomposition regularization method is as follows: Based on the previous two methods, we introduce a new regularizing filter function, which is defined as follows: where α > 0, 0 < μ ≤ kKk, 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.
(3) It is obvious that qðα, μÞ = 1, as α → 0 So, from Lemma 1, 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: 4 Computational and Mathematical Methods in Medicine Then, a modified regularization method can be constructed based on the regularization operator R α 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 We combine the singular system theory of the compact operator with the previous theorem and show the following result for the regularized solution.

Error Analysis of the Regularized Solution
Theorem 2. Let x + = ðK * KÞ v z ∈ RðK * KÞ v , z ∈ X, with kzk ≤ 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: where x + is the solution of Kx = y, x δ αðδÞ ≔ R α y δ is the approximation of x, and K * denotes the adjoint operator of K.

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.
with the right-hand side and the analysis solution given by yðtÞ = e t and xðtÞ = 1.
We discretize the integral equation by the compound trapezoidal formula and obtain the linear system.

Computational and Mathematical Methods in Medicine
where A ∈ R 100×100 and the error-free right-hand side b = ½yðt 0 Þ,⋯,yðt m Þ T ∈ R 100 . The associated contaminated vector b is given by the following formula.
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 equation AX =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 (kAX −bk 2 , kXk 2 ) 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].      Computational and Mathematical Methods in Medicine 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.

Computational and Mathematical Methods in Medicine
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.
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 Section 2, 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.

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 Figure 10.  For the Fredholm integral equation of the first kind (Equation (12)), set the right-hand side equal to 1, that is, jðtÞ = 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.

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 function jðtÞ = 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 The inversed results obtained by using two different regularization methods are shown in Figure 15 and the optimized release profiles are depicted in Figure 16. From these two figures, we can see that the performance of MRM is bet-ter than that of TRM. The mean square deviations of two methods are 0.3019 and 0.2164, respectively.

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, jðtÞ 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 jj δ ðtÞ − jðtÞj ≤ δ. It is our aim to solve the following perturbed equation.

Computational and Mathematical Methods in Medicine
We use the modified regularization method to solve the above perturbed equation for three different cases. Figure 17 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.

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 con-  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

13
Computational and Mathematical Methods in Medicine 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.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The author declares no competing interests.