Optimal controlling of anti-TGF-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β and anti-PDGF medicines for preventing pulmonary fibrosis

In the repair of injury, some transforming growth factor-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}βs (TGF-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}βs) and platelet-derived growth factors (PDGFs) bind to fibroblast receptors as ligands and cause the differentiation of fibroblasts into myofibroblasts. When the injury repair is repeated, the myofibroblasts proliferate excessively, forming fibrotic tissue. We goal to control myofibroblasts proliferation and apoptosis with anti-transforming growth factor-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β (anti-TGF-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β) and anti-platelet-derived growth factor (anti-PDGF) medicines. The novel optimal regulator control problem with two controls (medicines) is proposed to simulate how to the preventing pulmonary fibrosis. Idiopathic pulmonary fibrosis (IPF) consists of restoring a system of cells, protein, and tissue networks with injury and scar. Myofibroblasts proliferation back to its equilibrium position after it has been disturbed by abnormal repair. Thus, the optimal regulator control problem with a parabolic partial differential equation as a constraint, zero flux boundary, and given specific initial conditions, is considered. The myofibroblast diffusion equation stands as a governing dynamic system while the objective function is the summation of myofibroblast, anti-TGF-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β and anti-PDGF medicines for the fixed final time. Here, myofibroblast is a nonlinear state of time while anti-TGF-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β and anti-PDGF are two unknown control functions. In order to solve the corresponding problem a weighted Galerkin method is used. Firstly, we convert the myofibroblast diffusion equation to a system of ordinary differential equations using the Lagrangian interpolation polynomials defined at Gauss-Lobatto integration points. Secondly, by the calculus of variations, the optimal control problem is solved successfully using canonical Hamiltonian and extended Riccati equations. Numerical results are given, and the plots are depicted. Moreover, solutions to the problem in which there is no control are compared. Numerical results show that, over time, the myofibroblast increases and then remains constant when there is no control. In contrast, the current solution decreases and vanishes after 300 days by prescribing controller medicines for anti-TGF-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β and anti-PDGF. The optimal strategy proposed in this paper helps practitioners to reduce myofibroblasts by controlling both anti-TGF-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β and anti-PDGF medicines.


Optimal controlling of anti-TGF-β and anti-PDGF medicines for preventing pulmonary fibrosis Fatemeh Bahram Yazdroudi & Alaeddin Malek *
In the repair of injury, some transforming growth factor-β s (TGF-β s) and platelet-derived growth factors (PDGFs) bind to fibroblast receptors as ligands and cause the differentiation of fibroblasts into myofibroblasts.When the injury repair is repeated, the myofibroblasts proliferate excessively, forming fibrotic tissue.We goal to control myofibroblasts proliferation and apoptosis with anti-transforming growth factor-β (anti-TGF-β ) and anti-platelet-derived growth factor (anti-PDGF) medicines.The novel optimal regulator control problem with two controls (medicines) is proposed to simulate how to the preventing pulmonary fibrosis.Idiopathic pulmonary fibrosis (IPF) consists of restoring a system of cells, protein, and tissue networks with injury and scar.Myofibroblasts proliferation back to its equilibrium position after it has been disturbed by abnormal repair.Thus, the optimal regulator control problem with a parabolic partial differential equation as a constraint, zero flux boundary, and given specific initial conditions, is considered.The myofibroblast diffusion equation stands as a governing dynamic system while the objective function is the summation of myofibroblast, anti-TGF-β and anti-PDGF medicines for the fixed final time.Here, myofibroblast is a nonlinear state of time while anti-TGF-β and anti-PDGF are two unknown control functions.In order to solve the corresponding problem a weighted Galerkin method is used.Firstly, we convert the myofibroblast diffusion equation to a system of ordinary differential equations using the Lagrangian interpolation polynomials defined at Gauss-Lobatto integration points.Secondly, by the calculus of variations, the optimal control problem is solved successfully using canonical Hamiltonian and extended Riccati equations.Numerical results are given, and the plots are depicted.Moreover, solutions to the problem in which there is no control are compared.Numerical results show that, over time, the myofibroblast increases and then remains constant when there is no control.In contrast, the current solution decreases and vanishes after 300 days by prescribing controller medicines for anti-TGF-β and anti-PDGF.The optimal strategy proposed in this paper helps practitioners to reduce myofibroblasts by controlling both anti-TGF-β and anti-PDGF medicines.
After cell destruction, macrophages and other cells begin to produce inflammatory mediators (messenger molecules), including TGF-β and PDGF to transform fibroblasts to myofibroblasts 1,2 .Myofibroblast cells appear during wound repair.Myofibroblasts secrete large amounts of extracellular matrix (ECM).The activity of these cells causes wound closure after injury.TGF-β is a potent protein in enhancing collagen production by fibroblasts and myofibroblasts 3 .Moreover, PDGF proteins localized and sustained caused abnormal fibroblast proliferation and collagen production in IPF.There is PDGF protein production in macrophages and epithelial cells of patients but not in normal lung tissue 4 .
Hao first presented a mathematical model for sarcoidosis as a biomedical problem in 2014 5 , then by developing his model they used it for chronic pancreatitis 6 .A mathematical model of the interstitial fibrosis immune system is proposed by Hao et al. 7 They monitored the effectiveness of existing anti-fibrotic drugs or those undergoing clinical trials in renal fibrosis.M1-derived inflammatory macrophages and M2 anti-inflammatory alveolar macrophages were considered for pulmonary fibrosis problems 8 .Hao et al. used this model to evaluate the effect of other potential drugs aimed at preventing liver fibrosis in 2017 9 .
Optimal control theory is a significant branch of modern control.Deals with the best possible control strategy that minimizes a certain performance index.It can be used as a powerful tool to solve the optimal control problem of disease dynamics.Solving optimization problems subject to constraints given in terms of partial differential equations (PDEs) is one of the most challenging problems.However, in medical, industrial, and (i) A novel dynamic system is modeled during the time of inflammation and drug administration.(ii) A new hybrid optimal control problem with PDE constraint with two controls is applied.(iii) An optimal control problem for decreasing myofibroblast is proposed where both anti-TGF-β and anti-PDGF medicines are controlled by a novel technique.(iv) The spectral method is used for discretization.(v) The myofibroblast diffusion equation is converted to a system of ordinary differential equations using the Lagrangian interpolation polynomials defined at Gauss-Lobatto integration points.(vi) Canonical Hamiltonian and extended Riccati equations for two controls are used.(vii) the extended linear feedback is used to solve the related optimal control problem.(viii) A constant vector c appears in the related state space ordinary differential equation (see Eq. 57).
(ix) The optimal strategy proposes to control both anti-TGF-β and anti-PDGF medicines.
In the present paper, in Fig. 1, IPF is shown schematically.In Fig. 2, differentiate myofibroblast from fibroblast is shown.In Figs. 3 and 4, lung tissue with and without damage area is shown.The myofibroblast diffusion and homogenized equations are proposed.The functions, variables, and parameters are given in Table 1.Legendre polynomials and Gauss-Lobatto integration points for Galerkin spectral method are introduced.This optimal control problem is solved by first discretizing and then optimizing technique.Firstly, the myofibroblast diffusion partial differential equation (PDE) is converted into an algebraic system of ordinary differential equations (ODEs) by the Galerkin spectral method.Secondly, the optimal control problem using Pontryagins minimum principle is solved.Numerical results are given in Figs. 5, 6, 7, 8, 9,10, 11, 12 and Tables 2, 3 and 4. Finally, discussion and conclusions are presented.

Figure 1.
Part of schematic network for a cell, protein, and tissue in Idiopathic pulmonary fibrosis (IPF).When an injury occurs in an organ, the immune system secretes pro-inflammatory cytokines to suppress and respond to the injury.Inflammatory responses, if excessive, cause serious damage to the inflamed tissue 22 .

Mathematical configuration
Lung tissue simulation.Figure 1 shows part of the cell schematic network, proteins, and tissues in IPF.
We assume the lung tissue is a square with an edge size of 1 cm.The square is divided into small squares and is called T ε with an edge size of ε .A simple representation of the lung geometry with two dimensions of x and y is considered.In each small square, there is a concentric circle as alveolar air space ( A ε ).Alveolar tissue is shown between the squares and circles in Fig. 3.
We assume that ε is extremely small and close to zero then we face a homogenized alveolar tissue ( T ε /A ε ).In this case, we ignore the alveoli space in the square and call it R square.Therefore, lung tissue is just a square without alveolar space, as shown in Fig. 4(b).Tissue inflammation is a small square D in R (R = [0, 1] × [0, 1] ).For a mild case of IPF, we assume that D = 0.3 × 0.3 cm 222 .

Myofibroblast diffusion equation. According to the functions, variables, and parameters given in
Table 1 the myofibroblast diffusion equation for one-dimension in D q , q = x or y.The myofibroblast diffusion equation for one-dimension is as follows:

R
The homogenized lung tissue

D q
The homogenized damaged tissue in one-dimension, q = x or y A ε Alveolar air space extremely small Figure 5. Dynamical system for myofibroblast density (no medication involves).The system of ODEs (35)  are solved for mN (t) during 400 days.Eq. ( 11) is discretized using the Lagrangian interpolation polynomials of order 16, 24, and 32 defined at Gauss-Lobatto integration points for green, blue, and red colors respectively.It is observed that, over time, the myofibroblast density increase to a certain amount.After approximately 300 days it remains constant but never vanishes.This shows that myofibroblasts resist apoptosis in response to serious injury.This means that if one uses no medication the persistent myofibroblast repairing leads to tissue remodeling and fibrosis formation 1,2 .and η G (t) (two medicines are involved).Just Eq. ( 39) is solved for m(t) during 400 days.Eq. ( 36) is discretized using the Lagrangian spectrol method for N = 32 with different η T (t) and η G (t) for (η T (t), η G (t)) = (0, 0), (0.1, 0.1), (0.3, 0.3), (0.5, 0.5) .It is observed that, over time, the myofibroblast density increase to a certain value after approximately 300 days, and it remains constant.It is shown that although by increasing the rate of medications η T (t) and η G (t) , myofibroblast density decreases, however, the patient is not going to be cured even if the medication is increased.In general myofibroblast density has an increasing form up to some days it will stay constant after and it never vanishes thus it never removes.This process shows that if one changes the medication doses self-willed, we cannot expect that the myofibroblast density will disappear and apoptosis will happen.This experience encouraged the authors to go for simulating the whole problem as an optimal control problem.6), we first convert Eq. ( 38) to a linear form (39).For discretization, we use the Lagrangian interpolation polynomials of orders 16, 24, and 32 defined at Gauss-Lobatto integration points.Then, the optimal control problem Eqs.[(41, 42) with initial and boundary conditions ( 6)] for different N = 16, 24 , and 32 is solved.According to Table 3, the solution to the optimal control problem is more accurate for N = 32.
Vol:.( 1234567890) Table 3. Solutions for the optimal regulator control problem for Eqs.(41), 42) with initial and boundary conditions (6) for N = 16, 24 , and 32 at times t = 50, 100, 150 , and 250.Myofibroblast densities are computed for the center point of the region D ( x = y = 0.45).In this figure, we compare two methods (without and with control) for myofibroblast density.A dynamical system is solved by spectral method N = 32 (also see Fig. 5 red).It is observed that, over time, the myofibroblast density increase and then remains constant but never vanishes.The optimal regulator control problem is solved and optimal myofibroblast density is depicted by the blue line using the spectral method N = 32 (also see Fig. 8a).It is observed that by controlling both TGF-β and anti-PDGF the myofibroblast density decreases and vanishes after almost 300 days.
Figure 11.Absolute error dynamical system solutions.The absolute error for the solution of Eq. ( 35) with an approximation of the polynomial of degrees 32 and 24 is depicted by purple color.Moreover, the absolute error for the solution of Eq. ( 35) with an approximation of the polynomial of degrees 32 and 16 is depicted by black color.
Figure 12.Absolute error optimal regulator control problem solutions.The absolute error for the solution of Eqs.[(41),(42), and ( 6)] using the spectral method with an approximation of the polynomial of degrees 32 and 24 is depicted in magenta color.Moreover, the absolute error for the solution of Eqs.[(41),(42), and ( 6)] using the spectral method with an approximation of the polynomial of degrees 32 and 16 is depicted in olive color.
Vol:.( 1234567890) The first term fibroblast transforms into myofibroblast by TGF-β and PDGF 23,24 .In Eq. ( 1), for simplicity, we set thus, the dynamical equation myofibroblast diffusion with initial and boundary conditions ( 2) is as follows: The myofibroblast diffusion equation for two-dimension in damaged area D is as follows: where, ∂y 2 , and initial and boundary conditions are

We set
The homogenized myofibroblast diffusion equation.According to Jikov et al. 25 and Goel et al. 26 , the homogenized myofibroblast diffusion equation are and where γ = 127 343 , ∇2 = a ∂ 2 ∂x 2 , and a = 0.11 .We divide both sides of Eqs. ( 8) and ( 9) by γ .The homogenized myofi- broblast diffusion equation for one and two dimensions are as follows and Spectral method.This part of the paper includes the essential formulas for Legendre polynomials, and Legendre spectral method in one and two dimensions, together with the discretization technique [27][28][29][30] (2) ( m(x 0 , y 0 , t 0 ) = 8.5 × 10 −3 , (x 0 , y 0 , t 0 ) = (0.3, 0.3, 0) ∂m(t) ∂x = 0, x = 0.3, 0.6, ∂m(t) ∂y = 0, y = 0.3, 0.6.L k (ξ ) is even when k is even and odd when k is odd.If L k (ξ ) is normalized so that L k (1) = 1 .For each k, we get: where [k/2] denotes the integral part of k/2.L 0 (ξ ) = 1 and L 1 (ξ ) = ξ .For each pair of Legendre polynomials of degrees k and M, the following orthogonality property applies where δ kM is Kronecker's delta.The kth-degree Lobatto polynomial, L0 k , derives from the (k + 1)-degree Leg- endre polynomial, L k+1 , as Legendre and Lobatto polynomials can be calculated by the recursive relations 27 Legendre spectral method in one-dimension.Basis functions are the Lagrangian interpolation polynomials defined at Gauss-Lobatto integration points.We define the approximate the order N = 16, 24, 32 for myofibroblast m N (x, t) , as follows where, mj (x, t) is the discrete polynomial coefficient of m j (x, t) and φ j is the j th Lagrange polynomial of order N on the Gauss-Legendre-Lobatto (GLL) points {ξ j } N j=0 and in which L N and L ′ N are the Legendre polynomial of order N and its derivative, respectively.To convert the [−1, 1] to [a, b], we use the mapping function while inverse mapping yields.
For h = x b − x a , the stiffness 28 ( S q ), mass ( K q ) and constant coefficients ( C q ) matrices are as follows and q = x, y: Using the Gauss quadrature, we have 27 ( 12) where the GLL quadrature weights {w k } N k=1 are given in the following The mass matrix is diagonal when the nodal points are the same as the quadrature points since the Lagrange polynomials have cardinality properties 29,30 .
For q = x or y, we define H 1 (D q ) and H 1 0 (D q ) spaces as follows: For Eq. ( 10), proper approximation for m N (x, t) applies as a weighted Galerkin method.Find mN (q, t) ∈ H 1 0 (D q ) such that for all φ ∈ H 1 0 (D q ).
We apply the Green theory and get the weak form as follows From the boundary condition ( ∂m N (x, t) ∂x = 0), We substitute (18) in (30).Thus, mN (x, t) can be determined by solving the following ODE systems where the entries of the c m , S x , K x , and C x are defined in (3), ( 20), (21), and (22).
Legendre spectral method in two-dimensions.We assume that the domain considered is partitioned into the quadrilateral where is the reference square.The local approximating functions are the tensor product of the one-dimensional Legendre polynomials.The approximation of order N for the unknown function m N (x, y, t) in the reference square is as follows: The stiffness matrix S, the mass matrix K, and constant coefficients C matrices 28 , respectively, are defined as: where the entries of the S x and S y are defined in (20), the entries of the K x and K y are defined in (21), and the entries of the C x and C y are defined in (22).For (11), a suitable approximation for m N (x, y, t) applies as a weighted Galerkin method thus The homogenized myofibroblast diffusion with medicines dynamical system.After cell destruction and using medicines, the myofibroblast diffusion Eq. is (11) changed to where η T (t) is anti-TGF-β and η G (t) is anti-PDGF.It is clear that if η T (t) = η G (t) = 0 then Eq. ( 36) is equal to Eq. (11).
Optimal regulator control problem.The homogenized myofibroblast diffusion Eq. (36) using the initial and boundary conditions ( 6) stands as a governing dynamic system while the objective function is the summation of myofibroblast, anti-TGF-β and anti-PDGF medicines for the fixed final time.From now on, for simplicity, we use the following notions m(x, y, t) = m(t), T GF (x, y, t) = T GF (t) and G(x, y, t) = G(t).
where c is defined in Eq. ( 7).J(m(t), η T (t), η G (t), t) : R 2 × R 2 → R is the objective functional consists of two terms m(t) : R 2 → R (is the state function), (η parmeter, mfT is the activation rate of myofibroblast due to TGF-β , and mfG is activation rate of myofibroblast due to PDGF, d m is the death rate of myofibroblasts, K G is PDGF saturation for activation of myofibroblasts, K T GF is TGF-β saturation for apoptosis for alveolar tissue apoptosis, f is fibroblasts density, T GF (t) is the concentration of activated TGF-β at (x, y) position, and G is the concentration of activated PDGF at (x, y) position, t ∈ [0, 350] .
(the value of parameters are described in see Table 1).Note that Eq. ( 38) is a system that considers the time of inflammation and drug administration while Eq. ( 11) is not.Legendre spectral method is used to discretize Eq. ( 38).Thus we deal with the following ODEs (see Legendre spectral method).
where using ( 22), (34), and (35), we get A, B T , B G , and C b as follows Thus, the discrete optimal control problem is In the next, we apply Pontryagin's minimum principle 31,32 .
Pontryagin's minimum principle.The minimization of the performance index J will be done using Pontryagin's minimum principle.The extended Hamiltonian for (41) and ( 42) is where, (t) is the vector of the Lagrange multipliers.Define J by The first differential J with respect to the vectors m(t) , η T (t) and η G (t) are given by A necessary condition for the performance index J to a minimum is that δ J = 0 .Thus, the vectors m(t) , η T (t) and η G (t) must satisfy in the following equations ( 36)  2 and are positive definite, this is sufficient to guarantee that η * T (t) and η * G (t) causes H( m(t), η T (t), η G (t), (t), t) to be a local minimum.
We use the linear feedback form for finding the optimal control, that is, look for functions K T (t) and K g (t).
For the unknowns ρ T , ρ G , K T (t) and K g (t) as the feedback matrices, we assume that the vector of the Lagrange multiplier * (t) is linear in m * (t) , i.e. for the unknown p(t) and g(t) if we substitute Eq. (54) in Eq. (52), we have By comparing (53) and (55), we have By substitute Eq. (55) in Eq. (42), we have From differentiate (54) and using (46), we have Finally, if we subsitute (54) and (57) in Eq. (58) We get In Eq. ( 60), m * (t) and c are positive and not zero.Thus, the coefficient of m(t) and the second term must be equal to zero simultaneously.Therefore Eq. ( 60) reduces to two differential equations (developed Riccati equations) as follows (52)

Numerical results
Numerical results are done using Python programming software version 3.8, while the processor is AMD Ryzen 5 5500U.Numerical results are presented as follows: (i) Just the dynamical system solution (no medication involves) The dynamical system is solved by transforming the related PDE to a system of ODEs using Lagrangian interpolation polynomials of order 16, 24, and 32 defined at Gauss-Lobatto integration points.In Fig. 5, the dynamical system for myofibroblast density (11) with initial and boundary conditions ( 6) is solved.Myofibroblast against time is plotted and compared.(ii) Just dynamical system solutions with different constant scalar values for η T (t) and η G (t) (two medicines are involved) By keeping the functions η T (t) and η G (t) as constant scalar values (η T (t), η G (t)) = (0, 0), (0.1, 0.1), (0.3, 0.3), (0.5, 0.5) in the dynamical system (36) myofibroblast densi- ties with different dosages of medications are computed and in Fig. 6 solutions are plotted.As it is shown in a problem simulation just with the dynamical system without medicines (11) and for the dynamical system with medicines (36) one can not cure the patient in this way.The reason is that the fibroblast density never vanishes and therefore never removes.Thus, to decrease and vanish myofibroblast density, one needs to change the problem formulation (36) to an optimal control problem.To do this, in the next steps the optimal regulator control problem (41, 42) with initial and boundary conditions ( 6) is proposed.(iii) Optimal regulator control problem solution (two medications as controls are involved) In the existence of two controls (medications) using the First Discretize, Then Optimize technique the optimal control problem in (37, 38) is solved.In Fig. 7, the optimal control problem solutions for myofibroblast density with two controls [Eqs.(41, 42)), initial and boundary conditions ( 6)] with different N (16,24,32) is depicted.
In Table 3, the optimal control problem solutions by the Lagrangian interpolation polynomials with 16, 24, and 32 degrees of freedom at times t = 50, 100, 150, 250 are shown.From this Table, one can recognize that even for moderate degrees of freedom solutions are accurate to 7 decimal points.From Table 4 one can find that m * 32 (t) gives a more accurate solution.Since, it is observed that the solution for N = 32 is more efficient, thus from now on, we use N = 32 in all of the computations.In Fig. 8(a), solution for optimal control problem Eqs.(41, 42) and ( 6) for myofibroblast density in point x = y = 0.45 when two controls (anti-TGF-β and anti-PDGF) appear in the related dynamical system using the spectral method are depicted.Fig. 8(b) shows behavior of anti-TGF-β and anti-PDGF in the center of the regian D during 350 days.In Fig. 9, the behavior of myofibroblast density against position (x) and time is plotted.It is observed that the myofibroblast density in the existence of both η T (t) and η G (t) decreases and vanishes after 300 days.(iv) Comparison between the solutions without and with two controls Fig. 10 shows the comparison between solutions when the solution no medicine is used (just the dynamical system solution is plotted) with when two medications exist (optimal regular control problem solution is plotted).It is observed that modeling the problem just the dynamic system gives a solution to myofibroblast density that never vanishes and therefore apoptosis will happen while the solution for optimal regulator problem decreases and vanishes after almost 300 days.The numerical result in Fig. 10 this verifies that the authors guest realistic assumptions in changing the problem modeling from just a dynamic system to an optimal regulator control problem are correct.(v) Convergence of the dynamical system and convergence of the solution of optimal control problem Tables 2 and 4 shows the absolute error of the dynamical system and the optimal control solution.It is observed that the absolute error is decreased when both degrees of freedom and the number of iterations are increased.The absolute error approximation of myofibroblast density with the Lagrangian interpolation polynomials of orders 16, 24, and 32 are calculated.From Tables 2, 3 and 4 wan can see that even for N = 16 degrees of freedom solutions for both m(t) and m * (t) are accurate up to 7 decimal points.Numerical results are plotted in Figs.11 and 12.

Discussion and conclusions
IPF is a chronic progressive disease of unknown etiology with approximately 5000 new cases per year and 5-year survival.This rate of incidence and mortality is higher than many other cancers.Furthermore, there is no proven effective treatment for IPF 1,33 .In this article, the homogenized diffusion equation is used to describe the space of lung alveoli.For the first time, we have proposed a mathematical optimal control problem with two control for the treatment of IPF.Anti-TGF-β and anti-PDGF medicines in myofibroblast diffusion are controlled successfully.First, the dynamical system of myofibroblast diffusion is solved by Legendre spectral method, and it is shown that using the spectral approximation with 32 nodes can give the proper solution (see Fig. 5).Hence myofibroblasts resist apoptosis in response to serious injury, and persistent repairing leads to tissue remodeling and fibrosis formation.This means that without medication can not expect to cure the disease.Even without any specific strategy, if we give some medication to the patient, myofibroblast density will not vanish (see Fig. 6).In Fig. 6, we showed that with the change of η T (t) and η G (t) in the dynamic system, we can see the reduction (61) ṗ(t) + p(t)A − p(t)B T p(t)B T − p(t)B G p(t)B G + I + AP(t) = 0 ġ(t) − p(t)B T g(t)B T − p(t)B G g(t)B G + p(t) + Ag(t) = 0 of myofibroblasts density but it never vanishes, and it never the cure diseases.Some researchers 5,7-9 use just the dynamical system and claim that by adjusting the medication doses can cure diseases and control fibrosis.Here, we show that in this manner there is no way to force the myofibroblast density to vanish and it never removes.For this reason, the authors model the problem as an optimal regulator problem with two controls as anti-TGF-β and anti-PDGF medicines.Here, we improve one of the models of Bahram Yazdroudi and Malek 21 to achieve the goals presented in this paper.It is observed that the control functions (anti-TGF-β and anti-PDGF) decrease and then remain zero after almost 300 days.Hence, in repair tissue, fibroblasts vanish through apoptosis, and no formation of fibrosis tissue happens.The medicines are prescribed from a certain dose, then decrease and vanish over time.With this strategy, there is no need to prescribe medicines during the last days of the patient take cure duration and the disease will be cured.When comparing two strategies (without and with control) for myofibroblast density, we consider that when there is no control, the myofibroblast increase and then remain constant (failed apoptosis).When there is control, after almost 300 days of controlling both anti-TGF-β and anti-PDGF, the myofibroblast density decreases and then vanishes.For example, the medicine Pirfenidone has been identified as an anti-TGF-β 34 (a TGF-β inhibitor that blocks TGF-β activity) and Imatinib as an anti-PDGF therapy, (a PDGF inhibitor that blocks PDGF activity).Numerical results in this paper, corroborate the idea of vanishing myofibroblast density by medication.To control myofibroblast proliferation, myofibroblast contraction, and apoptosis 34,35 , prescription of both anti-TGF-β and anti-PDGF medicines including antibodies is proposed.By this strategy, apoptosis and reduced myofibroblasts density prevent the formation of collagen in ECM 34 .It is observed that with the passage of time and taking medication, the myofibroblast density becomes zero after about 300 days.The patient needs both medicines anti-TGF-β for about 155 days and anti-PDGF for about 270 days to treat fibrosis.Here, in objective functional for the optimal control problem, the dosage of treatment through the use of anti-TGF-β and anti-PDGF medicines are the same.In further strategy, the dosage of anti-TGF-β and anti-PDGF medicines can be assumed to be different.In the further works the authors are going to discuss the effect of anti-IL13, in TGF-β and fibroblast, and the effect of anti-T α , in M1 and M2.

Figure 2 .
Figure 2. Differentiation of fibroblasts into myofibroblasts.Fibroblasts are activated and differentiate into myofibroblasts in response to tissue injury to initiate repair.Myofibroblasts secrete large amounts of ECM for repairing and remodeling.In normal repair, myofibroblasts vanish through apoptosis.In response to serious injury, myofibroblasts resist apoptosis.Their persistence leads to tissue remodeling and the formation of fibrosis.

Figure 3 .
Figure 3. Lung geometry without damage area.Lung geometry consists of squares with smaller circles in the center that show the alveolar air space.

Figure 4 .
Figure 4. Lung geometry with damage area.In (a), the AA square is lung geometry consisting of lots of squares with small circles in the middle.Circles show the alveolar air space.In (b), the alveolar air space is not considered since circles are extremely small.A damaged area D is shown in the center of the square R. The boundary conditions for D have zero flux.

Figure 6 .
Figure 6.Dynamical system for myofibroblast density with different constant scalar values for η T (t) and η G (t) (two medicines are involved).Just Eq. (39) is solved for m(t) during 400 days.Eq. (36) is discretized using the Lagrangian spectrol method for N = 32 with different η T (t) and η G (t) for (η T (t), η G (t)) = (0, 0), (0.1, 0.1), (0.3, 0.3), (0.5, 0.5) .It is observed that, over time, the myofibroblast density increase to a certain value after approximately 300 days, and it remains constant.It is shown that although by increasing the rate of medications η T (t) and η G (t) , myofibroblast density decreases, however, the patient is not going to be cured even if the medication is increased.In general myofibroblast density has an increasing form up to some days it will stay constant after and it never vanishes thus it never removes.This process shows that if one changes the medication doses self-willed, we cannot expect that the myofibroblast density will disappear and apoptosis will happen.This experience encouraged the authors to go for simulating the whole problem as an optimal control problem.

Table 2 . 9 Figure 7 .
Figure 7. Novel optimal regulator control problem with different N.The optimal myofibroblast density with two controls is depicted.For solving Eqs.(41, 42) with initial and boundary conditions (6), we first convert Eq. (38) to a linear form (39).For discretization, we use the Lagrangian interpolation polynomials of orders 16, 24, and 32 defined at Gauss-Lobatto integration points.Then, the optimal control problem Eqs.[(41, 42) with initial and boundary conditions(6)] for different N = 16, 24 , and 32 is solved.According to Table3, the solution to the optimal control problem is more accurate for N = 32.

Figure 8 .
Figure8.Novel optimal regulator control problem.In (a), optimal myofibroblast density with two controls is depicted using the spectral method ( N = 32 ).For solving Eqs.(41, 42) with initial and boundary conditions (6), we first transform Eq. (38) to a linear form (39).For discretization, we use the Lagrangian interpolation polynomials of order 32 defined at Gauss-Lobatto integration points.Then, we solve optimal control problem Eqs.(41, 42) with initial and boundary conditions (6) using Pontryagin's minimum principle, Hamiltonian and developed Riccati equations.It is observed that when we control both anti-TGF-β and anti-PDGF the myofibroblast density vanishes after almost 300 days.This strategy can be applied by physicians when they prescribe anti-TGF-β and anti-PDGF medicines in almost 300 days.In (b), the optimal control functions η T (t) and η G (t) are depicted.It is observed that the control functions (anti-TGF-β and anti-PDGF) decrease and then remain zero.Hence, in repair tissue, myofibroblasts vanish through apoptosis, and no formation of fibrosis tissue happens.The medicines are prescribed in certain doses and decrease over time.With this strategy, there is no need to prescribe medicines during some last days of the patient take cure duration.

Figure 9 .
Figure 9. Novel optimal regulator control problem.Myofibroblast density against position (x) and time is plotted.It is observed that the myofibroblast density decreases and vanishes over time with controlling both η T (t) and η G (t).

Figure 10 .
Figure10.Comparison without and with control.In this figure, we compare two methods (without and with control) for myofibroblast density.A dynamical system is solved by spectral method N = 32 (also see Fig.5red).It is observed that, over time, the myofibroblast density increase and then remains constant but never vanishes.The optimal regulator control problem is solved and optimal myofibroblast density is depicted by the blue line using the spectral method N = 32 (also see Fig.8a).It is observed that by controlling both TGF-β and anti-PDGF the myofibroblast density decreases and vanishes after almost 300 days.
(t) ṗ(t) + p(t)A − p(t)B T p(t)B T − p(t)B G p(t)B G + I + Ap(t) + c − p(t)B T g(t)B T − p(t)B G g(t)B G + p(t)We solve (61) by the Euler approximation method and using (54).Firstly, we find p(t), g(t) and * (t) from (61 and 58).Secondly, the problem (41) and (42) with initial and boundary conditions (6) can be solve.