EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 4, 2024, 2651-2675 ISSN 1307-5543 – ejpam.com Published by New York Business Global Simultaneous Identification of the Parameters in the Mathematical Model of Brain Tumor Growth Dynamics Under Treatment Salih Tatar1,∗, Mohamed BenSalah2, Maryam Alamil1 1 College of Science and General Studies, Alfaisal University, Riyadh 11533, Saudi Arabia 2 Department of Computer Sciences, Higher Institute of Applied Science and Technology of Sousse, University of Sousse, Rue Tahar Ben Achour, Sousse 4003, Tunisia Abstract. This paper is devoted to an inverse problem for a nonlinear parabolic equation related to brain tumor dynamics. After reformulating the inverse problem as a minimization problem, we prove the existence and stability of the solution to the minimization problem. Based on the Fréchet differentiability of the objective (cost) functional, we develop an efficient iterative procedure for the numerical solution to the minimization problem. Numerical examples with noise-free and noisy data illustrate applicability and accuracy of the proposed method to some extent. 2020 Mathematics Subject Classifications: 35R30, 65N21, 35Q92, 65M32, 35K61 Key Words and Phrases: Brain Tumor, Inverse Problem, Optimization Problem, Fréchet differentiability Nomenclature g(c) Reaction term B Spatial domain of the brain ν The outward normal vector on the boundary φ Initial cell density D(x) Diffusion coefficient c Cell density at position x and time t div Divergence operator ∇ Gradient operator ∂B Boundary of B h Mesh step in x direction τ Time step in t direction T Final time ε Noise level clim Carrying capacity ρ Proliferation parameter dw Diffusion coefficient for white matter dg Diffusion coefficient for gray matter k(x, t) Treatment term α Spatial component of the treatment profile β Temporal component of the treatment profile Aα ×Aφ Set of admissible solutions J(α,φ) Objective functional ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v17i4.5386 Email addresses: statar@alfaisal.edu (S. Tatar), mohamed.bensalah@fsm.rnu.tn (M. BenSalah), malamil@alfaisal.edu (M. Alamil) https://www.ejpam.com 2651 Copyright: © 2024 The Author(s). (CC BY-NC 4.0) S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2652 ζkα Step size for α ζkφ Step size for φ Sk α Search direction for α Sk φ Search direction for φ ε̄ Tolerance parameter ωi Weight functions 1. Introduction In this paper, we consider an inverse problem for the following problem: ∂tc− div(D(x)∇c) + k(x, t) c− g(c) = 0 (x, t) ∈ B × (0, T ), D(x) ∂νc = 0 (x, t) ∈ B × (0, T ), c(x, 0) = φ(x) x ∈ B, (1) where B ⊂ Rd (d = 1,2,. . . ) be an open bounded domain with a sufficiently smooth boundary ∂B , T > 0 is a final time, c(x, t) is the tumor cell density at the time t and the spatial location x within the brain region B, k(x, t) c(x, t) is the treatment term describing the death of cells due to chemotherapy or radiation therapy, D(x) is the diffusion coeffi- cient of cells in the brain tissue modeling a random tumor cell movement as a diffusive flux proportional to the cell density gradient, div(D(x)∇c) captures random tumor cell movement, g(c) is the reaction term that describes the growth dynamics of the tumor and ν is the outward unit normal. The boundary condition D(x) ∂νc = 0 ensures tumor cells do not diffuse beyond the brain region. Additionally, the function φ(x) ∈ L2(B) represents the initial tumor cell density at t = 0. The problem (1) is a well-known model related to the growth of brain tumors under treatment The model (1) serves as a fundamental tool for understanding brain tumor dynamics and optimizing treatment strategies through comprehensive mathematical analysis. For instance, the authors work the problem (1) with a source term such that the model has logistic growth and accomplish full 3-dimensional simulations of the tumor in time on two types of imaging data, the 3d Shepp-Logan head phantom image and an MRI T1-weighted brain scan from the Internet Brain Segmentation Repository in [15]. In [1], the authors study a treatment parameter identification problem in (1). They derive a nonlinear con- jugate gradient method for the inverse problem and reconstruct the unknown parameter from additional information about the tumor taken at a fixed instance of time. In [26], the authors analyze the behavior of gliomas mechanistically, deriving their behavior from two fundamental properties, the net proliferation rate and the diffusion coefficient. They conclude that although these clearly determine the tumor’s growth as an expansion of its traveling wave front, other factors (such as age and Karnofsky performance score) combine to determine the patient’s prognosis, the survival time. The authors extend a mathematical model of gliomas based on proliferation and diffusion rates to incorporate the effects of augmented cell motility in white matter as compared to grey matter and they simulate model tumors on an anatomically accurate brain domain by using a de- tailed mapping of the white and grey matter in the brain developed for a MRI simulator, S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2653 in [24]. A simple spatio-temporal mathematical model, based on proliferation and diffu- sion, that incorporates the effects of radiotherapeutic and chemotherapeutic treatments has been studied in [19]. Moreover, they study the effects of different schedules of radiation therapy, including fractionated and hyperfractionated external beam radiotherapy, using a generalized linear quadratic model. In [23], the authors present a novel explicit triscale reaction-diffusion numerical model of glioblastoma multiforme tumor growth, they adopt a finite-difference time-domain method for the numerical solution and then a clinical sce- nario is addressed to demonstrate the workflow of a possible clinical validation procedure. Other researchers [20] present an extension of Swanson’s reaction-diffusion model in order to include the effects of radiation therapy using the classic linear-quadratic radiobiological model for radiation efficacy, along with an investigation of response to various therapy schedules and dose distributions on a virtual tumor. Some of the recent developments in mathematical modeling are reviewed in [25]. A model of untreated gliomas is repre- sented first followed by models of polyclonal gliomas that follows a chemotherapy or a surgical resection. Such reviewed models illustrate the evolution of mathematical models for glioma growth and invasion beginning in simple homogeneous tissue, with or without gross anatomical boundaries (skull and ventricles), extending to complex heterogeneous tissue, with varying proportions of grey and white matter in cerebral cortex (including the sulcal pattern), deep cerebral nuclei, brainstem and cerebellum. A patient-specific, biologically based mathematical model for glioma growth that quantifies response to XRT in individual patients in vivo is presented in [21]. This mathematical model uses net rates of proliferation and migration of malignant tumor cells to characterize the tumor’s growth and invasion along with the linear-quadratic model for the response to radiation therapy. In this paper, we study an inverse problem for the nonlinear parabolic problem (1). Inverse problems of finding coefficients, the initial distribution, sources or parameters are physically and practically very important. Recently, there has been a growing interest in inverse problems for nonlinear parabolic equations. For instance, the simultaneous reconstruction of the initial temperature and heat radiative coefficient in a heat conduc- tive system is studied in [28]. They prove the stability of the inverse problem and then the reconstruction process is done by Tikhonov regularization. In [3], the identification of the space-dependent reaction coefficient, the initial temperature and the source term from measured temperatures at two instants t1, t2 and at the final time T are investi- gated. In [4], an inverse problem of simultaneously identifying and reconstructing the space-dependent reaction coefficient and source term component from time-integral tem- perature measurements is investigated. After they prove the existence and uniqueness of the solution, the conjugate gradient method to find the numerical solution is devel- oped, and its convergence is proved from the Lipschitz continuity of these gradients. A novel inverse problem of reconstructing the unknown time-dependent source term enter- ing the fourth-order parabolic equation of thermal grooving by surface diffusion from a given integral observation is formulated in [5]. The author studies the determination of an unknown time-dependent source term in a Kuramoto-Sivashinsky equation from given additional integral-type measurement in [2]. In [22], some methods are developed for find- S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2654 ing the thermophysical parameters of a two-layer soil and a rational method of choosing the damping coefficient is also proposed, which provides an indicative rate of convergence of the approximate value of the functional to zero. Our paper can be considered as an addition to these studies. This article is organized as follows: In the next section, we formulate and analyze the direct and the inverse problems. In Section 3 , we present an iterative procedure based on the conjugate gradient method to solve the minimization problem. Some numerical examples with noise-free and noisy data are given to show the efficiency of the method in Section 4. The conclusions and possible directions on the problem are given in Section 5. 2. Formulations and analyses of the direct and the inverse problems In this section, we define the direct and the inverse problems. Throughout this paper, we define Q := B×(0, T ) and BT := ∂B×(0, T ). The diffusion coefficient D(x) is assumed to encompass two regions, white and grey matter, as in the papers [10, 24, 25], reflecting the spatial heterogeneity of brain tissue. Specifically, we define D(x) as: D(x) = { dw, if x ∈ white matter, dg, if x ∈ grey matter, where dg and dw denote the diffusion coefficients for grey and white matter, respectively, satisfying 0 < dg ≪ dw. For the growth function g(c), we adopt a logistic growth model to ensure that the growth rate diminishes as the cell density approaches its maximum capacity. Thus, g(c) is defined as: g(c) = ρ c ( 1− c clim ) , where ρ is the net proliferation rate and clim represents the carrying capacity. This for- mulation provides a more realistic representation of tumor growth dynamics compared to simple exponential models [7, 9, 24, 26, 27], which may result in unbounded growth over time. Definition 1. The weak solution c ∈ L2 ( 0, T ;H1(B) ) of the problem (1) is defined as the solution of the following problem:∫ T 0 ∫ B ( ∂tc(x, t)ψ+D(x)∇c(x, t) ·∇ψ+k(x, t) c(x, t)ψ ) dx dt = ∫ T 0 ∫ B g(c)ψ dx dt, (2) for all ψ ∈ L2 ( 0, T ;H1(B) ) . For given inputs D(x), k(x, t), φ(x) and T > 0, the problem (1) is called the direct problem. Like most direct problems of the mathematical physics, the problem (1) is well posed, see [1, 16] for details. S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2655 Here the inverse problem consists of simultaneously determining the treatment profile k(x, t) and the initial tumor cell density φ, by means of the time-integral measurements of tumor cell density (measured data) defined as follows: ϕ1 := ∫ T 0 ω1(t) c(x, t) dt and ϕ2 := ∫ T 0 ω2(t) c(x, t) dt, x ∈ B, (3) where ω1(t) and ω2(t) denote two distinct weight functions. We note that since time- integral measurements of tumor cell density do not typically reveal time-dependent profiles, we consider the function k(x, t) ∈ L∞(Q) in the following separable form: k(x, t) = α(x)β(t), (4) where α(x) ∈ L∞(B) represents the spatial component of the treatment profile k(x, t) and describes how the treatment efficacy varies across different spatial locations within the brain region. On the other hand, β(t) ∈ L∞(0, T ) is the temporal component of the treatment profile and captures how the treatment effectiveness evolves over time. In the inverse problem considered here, our focus shifts to scenarios where the spatial component α(x) is unknown, while the temporal component β(t) is known. By taking (4) into account in the problem (1), we have ∂tu− div(D(x)∇u) +α⋆(x)β(t)u− g(u) = 0 in Q, D(x) ∂νu = 0 on BT , u(x, 0) = φ⋆(x) in B. (5) The focus of our investigation in this article revolves around the simultaneous determi- nation of the response treatment parameter α⋆(x) and the initial cell density φ⋆(x) from two time-integral measurements of u defined by (3). We first define the sets of admissible solutions to analyze the direct problem (5). A set Aα ×Aφ defined by Aα := {α(x) : α ∈ L∞(Ω), 0 < α ⩽ α(x) ⩽ α, a.e x ∈ Ω}, Aφ := {φ(x) : φ ∈ L∞(Ω), |φ(x)| ⩽ φ, a.e x ∈ Ω}, is called the set of admissible solutions where α, α, and φ are predetermined positive constants. For each pair (α,φ) belonging to Aα×Aφ, we denote the solution to problem (5) by c(α,φ). In real-world scenarios, observed data ϕ1 and ϕ2 are given with some error including measurement errors, biological variability, spatial heterogeneity and temporal dynamics inherent in tumor growth and treatment response. To overcome these effects, we introduce perturbed measurements ϕε1 and ϕε2 as follows: ∥ϕεi − ϕi∥L∞(Ω) ⩽ ε, for i = 1, 2, (6) S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2656 where ε > 0 represents the magnitude of noise in the data. Consequently, the quasi- solution, see [14], of the inverse problem under consideration is obtained by minimizing the following least-squares objective functional: J(α,φ) := 1 2 ∥∥∥∫ T 0 ω1(t) c(α,φ) dt− ϕε1 ∥∥∥2 L2(Ω) + 1 2 ∥∥∥∫ T 0 ω2(t) c(α,φ) dt− ϕε2 ∥∥∥2 L2(Ω) . (7) Therefore, the inverse problem addressed in this article can be reformulated and modeled by the following minimization problem: Find (α⋆,φ⋆) ∈ Aα ×Aφ such that J(α⋆,φ⋆) ⩽ J(α,φ), ∀(α,φ) ∈ Aα ×Aφ. (8) It is worth noting that in this paper, we adopt the regularization before discretization strategy, as opposed to the alternative approach of discretization before regularization discussed in [3, 11]. Theorem 1. There exists at least one minimizer for the optimization problem (8). Proof. Since J(α,φ) is nonnegative, we know that inf J(α,φ) is finite over Aα ×Aφ. Thus, there exists a minimizing sequence {αn,φn}n⩾0 from Aα ×Aφ such that lim n→∞ J(αn,φn) = inf α∈Aα φ∈Aφ J(α,φ). This implies the boundedness of {αn,φn}n⩾0 in L∞(B) × L2(B) and therefore there is a subsequence (still denoted as {αn,φn}n⩾0) such that both {αn}n⩾0 and {φn}n⩾0 converge weakly to α⋆ ∈ L∞(B) and φ⋆ ∈ L2(B), respectively, since the sets Aα and Aφ are closed and convex. We shall prove first that (α⋆,φ⋆) is indeed a minimizer of (8). Since each pair (αn,φn) corresponds to a solution cn := c (αn,φn) to (5) with α = αn and φ = φn, by letting ψ = cn in (2) and using 1 2∂tc 2 = ctc, it follows immediately that the sequence {cn := c (αn,φn)} is also bounded in L2(0, T ;H1(B)). This indicates the existence of some c⋆ ∈ L2(0, T ;H1(B)) and a subsequence of {cn}, still denoted by {cn}, such that cn := c (αn,φn)⇀ c⋆ in L2(0, T ;H1(B)). (9) We claim c⋆ = c (α⋆,φ⋆). The continuity of g(c) with respect to c together with (9) imply ∂tc (α n,φn)⇀ ∂tc ⋆, ∇c (αn,φn)⇀ ∇c⋆ in L2(Q). Letting n→ ∞ in (2) with α = αn and φ = φn yields∫ T 0 ∫ B (∂tc ⋆ ψ +D(x)∇c⋆ · ∇ψ +α⋆(x)β(t) c⋆ ψ) dx dt = ∫ T 0 ∫ B g(c⋆)ψ dx dt, (10) S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2657 for all ψ ∈ L2 ( 0, T ;H1(B) ) . Next, we shall prove c⋆(·, 0) = φ⋆, which with (10) implies c⋆ = c(α⋆,φ⋆). In doing so, we take ϑ in C1[0, T ], with ϑ(0) ̸= 0 and ϑ(T ) = 0 and ζ ∈ L2(B). We employ the classical integration by parts formula to get∫ T 0 ∫ B ∂tc (α n,φn) ψ ζ dx dt = − ∫ B φn ψ(0) ζ dx− ∫ T 0 ∫ B c (αn,φn) ∂tψ ζ dx dt. (11) Letting n→ ∞ in (11), we conclude that∫ T 0 ∫ B ∂tc ⋆ ψ ζ dx dt = − ∫ B φ⋆ ψ(0) ζ dx− ∫ T 0 ∫ B c⋆ ∂tψ ζ dx dt. (12) On the other hand, we have:∫ T 0 ∫ B ∂tc ⋆ ψ ζ dx dt = − ∫ B c⋆(·, 0)ψ(0) ζ dx− ∫ T 0 ∫ B c⋆ ∂tψ ζ dx dt, (13) which, together with (12), implies that c⋆(·, 0) = φ⋆. Therefore, c⋆ := c(α⋆,φ⋆). By αn ⇀ α⋆ and φn ⇀ φ⋆ in L∞(B) ⊂ L2(B) and L2(B), respectively, we use the lower semi-continuity of the L2-norm to conclude that J (α⋆,φ⋆) = 1 2 ∥∥∥∫ T 0 ω1(t) c(α ⋆,φ⋆) dt− ϕε1 ∥∥∥2 L2(B) + 1 2 ∥∥∥∫ T 0 ω2(t) c(α ⋆,φ⋆) dt− ϕε2 ∥∥∥2 L2(B) ≤ lim inf n→∞ 1 2 ∥∥∥∫ T 0 ω1(t) c(α n,φn) dt− ϕε1 ∥∥∥2 L2(B) + 1 2 ∥∥∥∫ T 0 ω2(t) c(α n,φn) dt− ϕε2 ∥∥∥2 L2(B) ≤ lim inf n→∞ J (αn,φn) = inf α∈Aα φ∈Aφ J(α,φ). Hence, the proof is complete. After that, we prove the stability of (8) in the sense that the minimization problem (8) effectively stabilizes the considered inverse problem against perturbations in the measured data. For this purpose, let {ϕε1,ℓ, ϕε2,ℓ}∞ℓ=0 denote a sequence of measurements of {ϕε1, ϕε2} in L2(B)×L2(B). For each ℓ ∈ N, we denote the solution of the minimization problem by (αℓ,φℓ). Then, we have the following minimization problem: min α∈Aα φ∈Aφ Jℓ(α,φ), Jℓ(α,φ) := 1 2 2∑ k=1 ∥∥∥∫ T 0 ωk(t) c(α,φ) dt− ϕεk,ℓ ∥∥∥2 L2(B) , ℓ = 0, 1, . . . (14) The following theorem states the convergence of the sequence {αℓ,φℓ}ℓ⩾0 when ϕεk,ℓ −→ ϕεk in L2(B) as ℓ −→ ∞, for k = 1, 2. (15) S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2658 Theorem 2. If (15) holds, then the sequence {αℓ,φℓ}ℓ⩾0 of minimizers of the problems (14) converges weakly in L∞(B)×L2(B) to a minimizer, associated with the measured data {ϕε1, ϕε2}, of the optimization problem (8). Proof. The existence of each pair (αℓ,φℓ) is guaranteed by Theorem 1. By definition, we have: Jℓ (αℓ,φℓ) ≤ Jℓ(α,φ), ∀(α,φ) ∈ Aα ×Aφ, that implies the uniform boundedness of (αℓ,φℓ) in L ∞(B)×L2(B). Consequently, there exist α⋆ ∈ L∞(B), φ⋆ ∈ L2(B), and a subsequence (still denoted as {αℓ,φℓ}ℓ⩾0), such that: (αℓ,φℓ)⇀ (α⋆,φ⋆) in L∞(B)× L2(B) as ℓ→ ∞. It suffices to show that (α⋆,φ⋆) is indeed a minimizer of (8). Repeating the same argument as that in the proof of Theorem 1, we derive: c (αℓ,φℓ)⇀ c (α⋆,φ⋆) in L2 ( 0, T ;H1(B) ) as ℓ→ ∞, up to taking a further subsequence. Combining this convergence with (15), we obtain:∫ T 0 ωk(t) c (αℓ,φℓ)−ϕεk,ℓ dt ⇀ ∫ T 0 ωk(t) c (α ⋆,φ⋆)−ϕεk dt in L2(B) as ℓ→ ∞, for k = 1, 2. Therefore, for k = 1, 2, we have:∥∥∥∫ T 0 ωk(t) c(α ⋆,φ⋆) dt− ϕεk ∥∥∥2 L2(B) ≤ lim inf ℓ→∞ ∥∥∥∫ T 0 ωk(t) c(αℓ,φℓ) dt− ϕεk,ℓ ∥∥∥2 L2(B) . By using the lower semi-continuity of the L2-norm, we deduce that: J (α⋆,φ⋆) = 1 2 2∑ k=1 ∥∥∥∫ T 0 ωk(t) c(α ⋆,φ⋆) dt− ϕεk ∥∥∥2 L2(B) ≤ 1 2 lim inf ℓ→∞ 2∑ k=1 ∥∥∥∫ T 0 ωk(t) c(αℓ,φℓ) dt− ϕεk,ℓ ∥∥∥2 L2(B) ≤ 1 2 lim ℓ→∞ 2∑ k=1 ∥∥∥∫ T 0 ωk(t) c(α,φ) dt− ϕεk,ℓ ∥∥∥2 L2(B) = 1 2 2∑ k=1 ∥∥∥∫ T 0 ωk(t) c(α,φ) dt− ϕεk ∥∥∥2 L2(B) = J(α,φ), ∀(α,φ) ∈ Aα ×Aφ. Thus, the proof is complete. S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2659 3. Conjugate gradient method In this section, we develop an efficient iterative procedure based on the conjugate gradient method to solve the minimization problem (8). Like all effective iterative methods for optimization problems, our approach requires information about the derivatives of the objective function. We start by deriving the Fréchet derivatives J ′ α(α,φ) and J ′ φ(α,φ) of the objective functional J(α,φ) with respect to the unknown terms α and φ, respectively. Subsequently, we establish the conjugate gradient method, which utilizes these gradients to reconstruct the unknown coefficients simultaneously. We refer the readers to [13] for convergence of this algorithm. 3.1. Fréchet Derivatives Following the methodology in [3, 4, 17], we can reduce the computational costs for the Fréchet derivatives by introducing the adjoint system of (5), which represents a backward time diffusion equation: −∂tv − div(D(x)∇v) + (α(x)β(t)− g′c(α,φ)) v = 2∑ k=1 ωk(t) (∫ T 0 ωk(τ) c(α,φ) dτ − ϕεk ) in Q, D(x) ∂νv = 0 on BT , v(x, T ) = 0 in B, (16) where g′c(α,φ) represents the Fréchet derivative of g at c(α,φ). The system (16) is linear and it is a well-posed problem, see [6]. For simplicity, we denote c(α,φ) by c. Theorem 3. The objective functional J(α,φ) is Fréchet differentiable and the derivatives J ′ α(α,φ) and J ′ φ(α,φ) at (α,φ) ∈ Aα ×Aφ are given by:J ′ α(α,φ) = − ∫ T 0 β(t) c(x, t) v(x, t) dt, J ′ φ(α,φ) = v(x, 0), (17) where v is the solution to the adjoint problem (16). We need the following two sensitivity problems to prove theorem 3. Sensitivity problem associated with the treatment parameter. Assume c(x, t) is perturbed by ϵ ηα(x, t) and the therapy parameter α(x) is perturbed by ϵα0(x), where ϵ > 0 is a small number and α0 ∈ Aα. Substituting the perturbed cell density c(x, t)+ϵ ηα(x, t) and the perturbed therapy parameter α(x) + ϵα0(x) into the original problem (5) yields a perturbed problem. The resulting sensitivity problem is given by ∂tηα − div(D(x)∇ηα) +α(x)β(t) ηα − g′c ηα = −β(t)α0(x) c in Q, D ∂νηα = 0 on BT , ηα(x, 0) = 0 in Ω. (18) S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2660 Sensitivity problem associated with initial activation of the tissue. Similarly, assume c(x, t) is perturbed by ϵ ηφ(x, t) and the initial activation of the tissue, denoted as φ(x), is perturbed by ϵφ0(x), where ϵ > 0 is a small number and φ0 ∈ Aφ. Substituting the perturbed cell density c(x, t) + ϵ ηφ(x, t) and perturbed initial activation of the tissue φ(x) + ϵφ0(x) into the original problem (5) yields a perturbed problem. The resulting sensitivity problem is given by ∂tηφ − div(D(x)∇ηφ) +α(x)β(t) ηφ − g′c ηφ = 0 in Q, D ∂νηφ = 0 on BT , ηφ(x, 0) = φ0(x) in Ω. (19) Note that the problems (18) and (19) are linear, and their well-posedness is classical, see [6]. These sensitivity problems are important to determine the step sizes in the descent direction in the conjugate gradient method. Proof of theorem 3. Let ϵ > 0 and α0 ∈ Aα be given, then direct calculations yield J(α+ ϵα0,φ)− J(α,φ) = ϵ { 2∑ k=1 ∫ Q ηα(x, t)ωk(t) (∫ T 0 ωk(τ) c(x, τ) dτ − ϕεk ) dx dt } + ϵ2 2 { 2∑ k=1 (∥∥∥∫ T 0 ωk(t) ηα dt ∥∥∥2 2 )} . (20) Therefore, the directional derivative of J(α,φ) with respect to α in the direction α0(x) is given by J ′ α(α,φ) ·α0 = lim ϵ→0 J(α+ ϵα0,φ)− J(α,φ) ϵ = 2∑ k=1 ∫ Q ηα(x, t)ωk(t) (∫ T 0 ωk(τ) c(x, τ) dτ − ϕεk ) dx dt. Multiplying both sides of the first equation in (16) by ηα, integrating over Q, using the Green formula and the fact that D(x) ∂νv = 0 on BT , we obtain − ∫ Q ∂tv ηα dx dt+ ∫ Q D(x)∇v∇ηα dx dt+ ∫ Q (α(x)β(t)− g′c) v ηα dx dt = 2∑ k=1 ∫ Q ωk(t) ηα (∫ T 0 ωk(τ) c(x, τ) dτ − ϕεk ) dx dt. Applying the Green’s and the integration by parts formulas in the above identity and taking into account that ηα(x, 0) = 0 and v(x, T ) = 0 for all x ∈ Ω, we get S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2661 ∫ Q ( ∂tηα − div(D(x)∇ηα) +α(x)β(t) ηα − g′c ηα ) v dx dt = 3∑ k=1 ∫ Q ωk(t) ηα (∫ T 0 ωk(τ) c(x, τ) dτ − ϕεk ) dx dt. The first equation in (18) leads to 2∑ k=1 ∫ Q ωk(t) ηα (∫ T 0 ωk(τ) c(x, τ) dτ − ϕεk ) dx dt = − ∫ Q β(t)α0(x) c(x, t) v(x, t) dx dt. Therefore, the directional derivative of J(α,φ) with respect to α in the direction α0(x) is given by J ′ α(α,φ) ·α0 = − ∫ Q β(t)α0(x) c(x, t) v(x, t) dx dt. (21) Similarly, by using the weak formulations of the adjoint problem (16) and the sensitivity problem (19), we can derive J ′ φ(α,φ) ·φ0 = ∫ Ω φ0(x) v(x, 0) dx. The proof is complete. 3.2. Iterative algorithm In this subsection, we employ a conjugate gradient method to approximate the mini- mizer of the functional J(α,φ) defined by (7). Let (αk,φk) be the kth approximation of the solution (α,φ). The iterations are given by the following equations: αk+1 = αk + ζkαS k α and φk+1 = φk + ζkφS k φ, k = 0, 1, 2, · · · . (22) In (22), the index k denotes the number of iterations, α0(x) and φ0(x) are the initial guesses for α(x) and φ(x), respectively. The terms ζkα and ζkφ are the step sizes, and Sk α and Sk φ represent the search directions which are determined by the following equations: Sk α = { −J ′0 α , k = 0, −J ′k α + ϑkα S k−1 α , k ≥ 1, and Sk φ = { −J ′0 φ , k = 0, −J ′k φ + ϑkφ S k−1 φ , k ≥ 1, (23) where J ′k α = J ′ α(α k,φk) and J ′k φ = J ′ φ(α k,φk) denote the derivatives of J(α,φ) with respect to α and φ, respectively, at point (αk,φk). The terms ϑkα and ϑkφ are the conjugate coefficients obtained using the Fletcher-Reeves formula [12] ϑkα = ∥∥J ′k α ∥∥ 2∥∥∥J ′k−1 α ∥∥∥ 2 and ϑkφ = ∥∥J ′k φ ∥∥ 2∥∥∥J ′k−1 φ ∥∥∥ 2 , k = 1, 2, · · · . (24) S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2662 Building upon the preceding discussions, it is noteworthy that all parameters are explicitly expressed, except the search step sizes ζkα and ζkφ. The determination of these parameters is accomplished through the utilization of the line search method. Specifically, the values for ζkα and ζkφ are obtained by minimizing the following functional: J(αk+1,φk+1) := 1 2 2∑ i=1 (∥∥∥∫ T 0 ωi(t)u(α k+1,φk+1) dt− ϕεi ∥∥∥2 2 ) . Since the search step sizes ζkα and ζkφ are implicit in the expression of J(αk+1, φk+1), we adopt a similar approach to [3, 4] and linearize it such that these parameters become explicit in the new formulation. Hence c(αk+1,φk+1) becomes as follows: c(αk+1,φk+1) :≈ c(αk,φk) + ζkαη k α + ζkφη k φ, (25) where ηkα and ηkφ are obtained by solving the sensitivity problems (18) and (19) with (α,φ) = (αk,φk) and (α0,φ0) = (Sk α, S k φ). We now define cki = ∫ T 0 ωi c k dt, ηkα,i = ∫ T 0 ωi η k α dt and ηkφ,i = ∫ T 0 ωi η k φ dt, for i = 1, 2, with ck = c(αk,φk). Then, from (7) and (25), we have J(αk+1,φk+1) = 1 2 2∑ i=1 ∫ Ω ( cki + ζkα η k α,i + ζkφ η k φ,i − ϕεi )2 dx. (26) It follows ∂J(αk+1,φk+1) ∂ζkα := R1ζ k α +R2ζ k φ − Y1 and ∂J(αk+1,φk+1) ∂ζkγ := R2ζ k α +R3ζ k φ − Y2 , where R1 = 2∑ i=1 ∥ηkα,i∥22, R2 = 2∑ i=1 ∫ B ηkα,i η k φ,i dx, R3 = 2∑ i=1 ∥ηkφ,i∥22, Y1 = 2∑ i=1 ∫ B ηkα,i(ϕ ε i − cki ) dx, and Y2 = 2∑ i=1 ∫ B ηkφ,i(ϕ ε i − cki ) dx. By setting ∂J(αk+1,φk+1) ∂ζkα = ∂J(αk+1,φk+1) ∂ζkφ = 0, one can deduce that the search step sizes ζkα and ζkφ satisfy the following linear system{ R1ζ k α +R2ζ k φ = Y1, R2ζ k α +R3ζ k φ = Y2. S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2663 Therefore, direct calculations yield ζkα = R3Y1 −R2Y2 R3R1 −R2 2 and ζkα = R1Y2 −R2Y1 R3R1 −R2 2 . (27) In the context of iterative algorithms, it is widely recognized that determining a suitable stopping criterion is crucial. Addressing this concern, the present study employs the discrepancy principle to conclude the iteration process. According to this principle, the iteration procedure if the following condition is fulfilled: J ( αk,φk ) ≤ ε̄, (28) where ε̄ is a small positive value, such as ε̄ = 10−6 for exact measurements. In cases where the measurements exhibit noise, ε̄ = 1 2 2∑ i=1 ∥ϕεi − ϕi∥22 . Building upon (6), it is established that ε̄ ⩽ ε2. Consequently, the primary steps of our reconstruction approach are summarized in the following algorithm: Algorithm 1 (H).Step 1. Set k = 0 and choose initial guesses α0 and φ0 for the unknown coefficients α and φ, respectively. Step 2. Solve the direct problem (5) with (α,φ) = (αk,φk) to get ck = c(αk,φk). Step 3. Solve the adjoint problem (16) and evaluate the gradients J ′,k α = J ′ α(α k,φk) and J ′,k φ = J ′ φ(α k,φk) given in theorem 3. Step 4. Calculate the conjugate coefficients ϑkα and ϑkφ by (24) and the directions Sk α and Sk φ by (23). Step 5. Calculate the step sizes ζkα and ζkφ by (27). Step 6. Update αk+1 and φk+1 by (22). Step 7 If the condition (28) is satisfied, then go to Step 8. Otherwise set k = k+1 and go to Step 2. Step 8 End. S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2664 4. Numerical experiments We set T = 1 and B = (0, 1) in all numerical examples. We rely on the finite difference method as our primary tool to deal with not only the direct problem but also its associated sensitive and adjoint variants. Setting a structured approach, we designate h = 0.02 as the step size and τ = 0.005 as the time increment, ensuring a balanced exploration of the problem domain both spatially and temporally. Moreover, to handle integrals within our computations, we employ the trapezoidal rule, a widely recognized numerical technique that provides reliable approximations. All numerical experiments are carried out in MAT- LAB (version R2017a). In the context of noisy observations, we enhance the fidelity of our simulations by incorporating Gaussian additive noise to ϕε1 and ϕε2, as detailed in (6). This noise, with amplitude σ, is dynamically adjusted based on a specified percentage p representing the noise level in the data. Utilizing the term ”random(1)”, we introduce variability represen- tative of real-world conditions by generating random values from a Gaussian distribution with a mean of zero and a standard deviation of unity. Specifically, the noisy observations ϕεi are generated according to: ϕεi = ϕi + σ × random(1), i = 1, 2, where σ is calculated as σ = p 100 × maxx∈Ω{|ϕ1(x)|, |ϕ2(x)|}. This approach ensures that our computational model captures the inherent uncertainties present in real-world measurements, enhancing the robustness of our simulations. The accuracy errors Eα(k) and Eφ(k), as functions of the iteration number k, for the treatment parameter α(x) and the initial cell density φ(x), respectively, are defined as Eα(k) := ∥αk −α∥2, (29) Eφ(k) := ∥φk −φ∥2, (30) where αk and φk are the CGM iterates, and α and φ are the analytical expressions of the treatment parameter and the initial cell density. Throughout this section, we maintain a consistent temporal component β(t) within the treatment profile, specified as exp (−t). Furthermore, our initial guesses for the treatment parameter α0 and initial cell density φ0 are both set to 1, establishing a standardized starting point for our analysis. Expanding upon our problem configuration, the parameters governing the logistic reaction function g(c) and the diffusion term are listed in Table 1. Parameter Value References Growth rate ρ 0.0012 Cock et al. [8] Diffusion coefficient dg 0.0013 Tracqui et al. [27] Diffusion coefficient dω 0.0065 Swanson et al. [24] Carrying capacity clim 0.625 Noviantri et al. [18] Table 1: Numerical values of physical coefficients. S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2665 Example 1. In this example, we aim to evaluate the performance of our algorithm in identifying the treatment parameter α and initial cell density φ within a controlled setting. These two functions are assumed to exhibit smooth behavior, reflecting typical characteristics observed in many practical applications. By selecting smooth functions for α and φ, we aim to establish a baseline assessment of the algorithm’s capability to accurately reconstruct these essential parameters under ideal conditions. This test serves as a crucial step in validating the effectiveness and reliability of our computational approach before extending it to more complex scenarios. In doing so, we choose the weight functions as ω1(t) = 1 and ω2(t) = e−t, tailoring our algorithm to specific characteristics of the problem domain. With these choices, we proceed to test our algorithm’s ability to reconstruct a treatment parameter α and initial cell density φ defined as α(x) = 2x+ sin(πx) + 1 and φ(x) = cos(4πx) + 1 2 . In Figure 1, we show the behavior of the objective functional J(αk,φk), as defined in (7), with respect to the iteration number k. This functional characterizes the simultaneous recovery of the two unknown parameters α(x) and φ(x) in the noise free case (p = 0), as well as with noise levels p = 1 and p = 2. From this figure, it can be seen that the objective functional (7) is monotonic decreasing and convergent. Additionally, we observe rapid convergence of the functional to a small positive value, underscoring the algorithm’s efficiency in minimizing the objective and thereby refining the parameter estimates. Figure 1: The objective functional J(αk,φk), for different noise levels p ∈ {0, 1, 2}, in Example 1. Based on the discrepancy principle (28), Table 2 presents a comprehensive overview of key metrics for each chosen noise level p ∈ {0, 1, 2}. Specifically, the table lists the stopping iteration numbers k⋆, the threshold ε, the errors (29) associated with the unknowns α(x) S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2666 and φ(x), as well as the norms of the Fréchet gradients ∥J ′,k⋆ α ∥2 and ∥J ′,k⋆ φ ∥2. ε̄ k⋆ Eα(k ⋆) Eφ(k ⋆) ∥J ′,k⋆ α ∥2 ∥J ′,k⋆ φ ∥2 p = 0 1.00e− 06 161 4.41e− 02 3.11e− 02 2.09e− 04 3.34e− 04 p = 1 8.54e− 05 66 1.21e− 01 1.17e− 01 3.04e− 04 3.63e− 04 p = 2 3.41e− 04 43 2.21e− 01 1.56e− 01 1.24e− 03 1.45e− 03 Table 2: Choices of noise levels p and their respective threshold values, the corresponding stopping iteration numbers k⋆, L2-errors, and norms of the Fréchet gradients in Example 1. From Table 2, it is evident that the numerical solutions are reasonably accurate for both α(x) and φ(x). Furthermore, the norms of the Fréchet gradients indicate that the Conjugate Gradient Method (CGM), with iterations stopped by the discrepancy principle (28), serves as a semi-convergent regularization method. In Figure 2, we illustrate the numerical solutions for the treatment parameter α(x) and the initial cell density φ(x) at the stopping iteration numbers k⋆ determined from Table 2. These solutions are shown for noise levels p in {0, 1, 2}. From Figure 2, it is evident that for noise-free data (p = 0), the analytical and numerical solutions overlap, making them indistinguishable graphically. However, as the level of noise increases, the accuracy of the solutions decreases. Overall, it can be observed that stable and accurate solutions are obtained for both coefficients α(x) and φ(x). (a) The treatment parameter α(x) (b) The initial cell density φ(x) Figure 2: The numerical results for Example 1 for p ∈ {0, 1, 2}. Left: the treatment parameter α(x). Right: the initial cell density φ(x). Example 2. In this example, we test our algorithm to recover a smooth treatment parameter and a non-smooth initial cell density, driven by considerations in brain tumor growth dynamics. Typically, treatment parameters in such scenarios exhibit smooth varia- tions, while initial cell densities may feature sharp changes or irregularities. By addressing these characteristics, we aim to enhance the applicability of our algorithm in real-world S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2667 medical contexts, where accurate characterization of tumor growth dynamics is crucial for effective treatment planning and monitoring. Moreover, investigating the algorithm’s performance in recovering both smooth and non-smooth parameters provides valuable in- sights into its robustness and versatility across diverse tumor growth patterns. To demonstrate the effectiveness of our proposed iterative method in this scenario, we define the exact treatment parameter α(x) and initial cell density φ(x) to be reconstructed as follows: α(x) = 1 + sin(4πx) 2 and φ(x) =  0, if x ∈ [0, 0.2] ∪ [0.8, 0.1], 2 sin(πx)− 1, if x ∈ [0.2, 0.5], −10x 3 + 8 3 , Otherwise. The measured data ϕεi , for i = 1, 2, are generated using the weight functions ω1(t) = 1+ t2 and ω2(t) = exp(−2t). Figure 3 illustrates the objective functional J(αk,φk) given by (7) for the simultaneous reconstruction of the unknown coefficients with noise levels p ∈ {0, 1, 2}. From this figure, it is evident that the objective function rapidly decreases monotonically and thus it is convergent. Figure 3: The objective functional J(αk,φk), for different noise levels p ∈ {0, 1, 2}, in Example 2. The choices of p in the test and the corresponding numerical performances are listed in table 3. S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2668 ε̄ k⋆ Eα(k ⋆) Eφ(k ⋆) ∥J ′,k⋆ α ∥2 ∥J ′,k⋆ φ ∥2 p = 0 1.00e− 06 167 1.14e− 02 1.33e− 01 3.43e− 04 2.58e− 05 p = 1 3.18e− 02 17 9.05e− 02 2.16e− 01 1.17e− 01 1.01e− 01 p = 2 1.27e− 01 9 1.06e− 01 3.07e− 01 3.30e− 01 1.22e− 01 Table 3: Choices of noise levels p alongside their respective threshold values, the corre- sponding stopping iteration numbers k⋆, L2-errors, and norms of the Fréchet gradients in Example 2. In comparison with Example 1, Table 3 highlights the impact of the non-smoothness of the initial cell density on the accuracy of our method. The corresponding numerical solutions for the treatment parameter α(x) and the initial cell density φ(x) at the stopping iteration numbers obtained based on the discrepancy principle (28) are illustrated in Figure 4. (a) The treatment parameter α(x) (b) The initial cell density φ(x) Figure 4: The numerical results for Example 2 for different noise levels p ∈ {0, 1, 2}. Left: the treatment parameter α(x). Right: the initial cell density φ(x). Similar to the previous test, Figure 4 demonstrates that the efficiency of our algorithm diminishes with increasing noise level p. Specifically, the numerical results exhibit high accuracy when no noise is present (0% noise added in the measured data ϕi, for i = 1, 2). It is noteworthy that the reconstruction of smooth functions displays higher accuracy compared to non-smooth cases. However, stable and accurate solutions are attained for both coefficients overall. Example 3. In this example, we test our algorithm to identify a non-smooth treat- ment parameter and a smooth initial cell density. This experimental setup stems from insights gained in brain tumor dynamics, where treatments such as radiation therapy or chemotherapy often induce localized variations in treatment effectiveness, resulting in non-smooth patterns in the treatment parameter. Conversely, the initial distribution of healthy and tumor cells within the brain tissue tends to exhibit a more uniform or grad- S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2669 ual pattern, reflecting the inherent composition and structure of the brain. By exploring this contrast in parameter characteristics, we aim to further enhance the adaptability and effectiveness of our algorithm in diverse tumor growth scenarios. For this purpose, we set ω1(t) = log(1+ t2) and ω2(t) = 2, and evaluate the algorithm’s performance in recovering: α(x) = max(1, 1− (1− x2) cos(5πx)) and φ(x) = cos(πx)− x2 + 2x 2 . The objective functional J(αk,φk) defined by (7) is shown in Figure 5. It is clear that the objective function decreases steadily and monotonically with the iteration numbers k and thus it is convergent. Figure 5: The objective functional J(αk,φk), for different noise levels p ∈ {0, 1, 2}, in Example 3. Based on the discrepancy principle (28) for each chosen noise level p ∈ {0, 1, 2}, we present the choices of p in this test along with the corresponding numerical performances in Table 4. ε̄ k⋆ Eα(k ⋆) Eφ(k ⋆) ∥J ′,k⋆ α ∥2 ∥J ′,k⋆ φ ∥2 p = 0 1.00e− 06 162 6.75e− 02 2.46e− 02 6.36e− 05 1.01e− 04 p = 1 1.83e− 05 116 7.92e− 02 5.71e− 02 6.62e− 04 6.94e− 04 p = 2 7.35e− 05 88 1.08e− 01 7.27e− 02 1.06e− 03 5.77e− 03 Table 4: Choices of noise levels p and their respective threshold values, the corresponding stopping iteration numbers k⋆, L2-errors, and norms of the Fréchet gradients in Example 3. Similar to the previous scenarios, we conclude this test by illustrating the comparisons of the recovered solutions with the exact ones in Figure 6. S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2670 (a) The treatment parameter α(x) (b) The initial cell density φ(x) Figure 6: The numerical results for Example 3 for different noise levels p ∈ {0, 1, 2}. Left: the treatment parameter α(x). Right: the initial cell density φ(x). As it is expected from the previous test, Figure 6 confirms that the efficiency of our algorithm is sensitive to both the noise level and the smoothness of the coefficient being reconstructed. Example 4. In brain tumor growth, various factors such as genetic mutations, cellular interactions, and microenvironmental conditions can contribute to the formation of com- plex spatial patterns in both the treatment response and the distribution of tumor cells. These patterns may include regions of high and low treatment effectiveness, as well as areas of dense and sparse tumor cell populations, resulting in non-smooth variations in these coefficients. According to this observation, in this test, we focus on the identification of the treatment parameter and the initial cell density in more complex scenarios, where both coefficients to be reconstructed are represented by non - smooth functions. In doing so, we choose the weight functions as ω1(t) = 1 + t + t2 and ω2(t) = t3, and we test our iterative procedure to identify the following two coefficients: α(x) = 2− |1− 2x| and φ(x) =  0, if x ∈ [0, 0.2] ∪ [0.8, 0.1], 10x 3 − 2 3 , if x ∈ [0.2, 0.5], −10x 3 + 8 3 , if x ∈ [0.5, 0.8]. Following the same steps as in the previous tests, we start this analysis by illustrating the evolution of the objective functional J(αk,φk) in Figure 7. This depicts the simultaneous recovery of the treatment parameter α(x) and the initial cell density φ(x), considering noise levels p = {0, 1, 2}. Notably, the objective function exhibits a monotonically de- creasing trend with k, converging rapidly to a small positive value, as expected. S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2671 Figure 7: The objective functional J(αk,φk), for different noise levels p ∈ {0, 1, 2}, in Example 4. The choices of noise levels p along with their corresponding numerical performances are detailed in Table 5. ε̄ k⋆ Eα(k ⋆) Eφ(k ⋆) ∥J ′,k⋆ α ∥2 ∥J ′,k⋆ φ ∥2 p = 0 1.00e− 06 207 6.23e− 02 1.52e− 02 6.21e− 04 1.61e− 04 p = 1 1.15e− 03 20 1.11e− 01 4.90e− 02 1.01e− 03 7.27e− 03 p = 2 4.60e− 03 11 1.58e− 01 7.38e− 02 1.09e− 03 6.52e− 03 Table 5: Choices of noise levels p alongside their respective threshold values, the corre- sponding stopping iteration numbers k⋆, L2-errors, and norms of the Fréchet gradients in Example 4. Utilizing the stopping iteration numbers k⋆, obtained according to the discrepancy principle (28), we present the comparisons between the exact solutions and the recovered ones for each chosen noise level p in Figure 8. S. Tatar, M. BenSalah, M. Alamil / Eur. J. Pure Appl. Math, 17 (4) (2024), 2651-2675 2672 (a) The treatment parameter α(x) (b) The initial cell density φ(x) Figure 8: The numerical results for Example 4 for different noise levels p ∈ {0, 1, 2}. Left: the treatment parameter α(x). Right: the initial cell density φ(x). In comparison with all previous tests, it is evident from Figure 8 that the accuracy and efficiency are reduced in this scenario due to the non-smoothness of the functions to be reconstructed. However, overall, the obtained results are acceptable. 5. Conclusions In this paper, we addressed an inverse problem involving the simultaneous determina- tion of the treatment response parameter and the initial tumor cell density in a nonlinear parabolic problem related to brain tumor dynamics. By reformulating the inverse problem as a minimization problem, we established the existence and initial stability of the solution. We proved the Fréchet differentiability of the objective (cost) functional and presented ex- plicit formulas for the derivatives by using the solution to related adjoint problem. Based on the Fréchet differentiability of the objective (cost) functional, we developed an effi- cient iterative procedure using the conjugate gradient method combined with Morozov’s discrepancy principle to solve the variational problem. Numerical examples, with both noise-free and noisy data, demonstrated the applicability and accuracy of the proposed method. On the other hand, several mathematical issues of high interest have not been discussed in this paper. The stability problem is one of them. The full stability issue is, however, to the best of our knowledge, still an open problem which deserves further inves- tigation. Additionally, the uniqueness analysis of the considered inverse problem remains an area for future research. These aspects will be the focus of forthcoming studies to enhance the robustness and reliability of the proposed methods in practical applications. Acknowledgements The authors thank to anonymous referees and editor for their valuable comments and suggestions that have greatly improved the presentation of the paper. REFERENCES 2673 References [1] G. Baravdish, B. T. Johansson, O. Svensson, and W. Ssebunjo. Identifying a response parameter in a model of brain tumor evolution under therapy. IMA Journal of Applied Mathematics, 88(2):378–404, 2023. [2] K. Cao. Identification of the time-dependent source term in a kuramoto–sivashinsky equation. Journal of Inverse and Ill-posed Problems, 32(3):361–387, 2024. [3] K. Cao and D. Lesnic. Simultaneous reconstruction of the spatially-distributed reac- tion coefficient, initial temperature and heat source from temperature measurements at different times. Computers & Mathematics with Applications, 78(10):3237–3249, 2019. [4] K. Cao and D. Lesnic. Simultaneous identification and reconstruction of the space- dependent reaction coefficient and source term. Journal of Inverse and Ill-posed Problems, 29(6):867–894, 2021. [5] K. Cao and D. Lesnic. Reconstruction of the time-dependent source in thermal grooving by surface diffusion. Journal of Computational and Applied Mathematics, 444:115789, 2024. [6] M. Choulli. Une introduction aux problèmes inverses elliptiques et paraboliques. Springer, 2009. [7] O. Clatz, M. Sermesant, P. Y. Bondiau, H. Delingette, S. K. Warfield, G. Malandain, and N. Ayache. Realistic simulation of the 3-d growth of brain tumors in mr images coupling diffusion with biomechanical deformation. IEEE transactions on medical imaging, 24(10):1334–1346, 2005. [8] J. Cook, D. E. Woodward, P. Tracqui, and J. D. Murray. Resection of gliomas and life expectancy. Journal of Neuro-Oncology, 24(4):131–135, 1995. [9] G. C. Cruywagen, D. E. Woodward, P. Tracqui, G. T. Bartoo, J. D. Murray, and EC Alvord. The modelling of diffusive tumours. Journal of Biological Systems, 3(4):937–945, 1995. [10] M. Dolgushin, V. Kornienko, and I. Pronin. Brain metastases: Advanced neuroimag- ing. Springer, 2018. [11] H. W. Engl and C. W. Groetsch. Projection-regularization methods for linear operator equations of the first kind. Australian National University, Mathematical Sciences Institute, 17:17–31, 1988. [12] R. Fletcher and C. M. Reeves. Function minimization by conjugate gradients. The computer journal, 7(2):149–154, 1964. REFERENCES 2674 [13] J. C. Gilbert and J. Nocedal. Global convergence properties of conjugate gradient methods for optimization. SIAM Journal on Optimization, 2:21–42, 1992. [14] V. K. Ivanov, V. V. Vasin, and V. P. Tanana. Theory of linear ill-posed problems and its applications. 2022. [15] R. Jaroudi, F. Aström, B. T. Johansson, and G. Baravdish. Numerical simulations in 3-dimensions of reaction–diffusion models for brain tumor growth. International Journal of Computer Mathematics, 2019. [16] B. Kaltenbacher and W. Rundell. On the simultaneous recovery of the conductiv- ity and the nonlinear reaction term in a parabolic equation. Inverse Problems and Imaging, 14(5):939–966, 2020. [17] Y. Liu, D. Jiang, and M. Yamamoto. Inverse source problem for a double hyperbolic equation describing the three-dimensional time cone model. SIAM Journal on Applied Mathematics, 75(6):2610–2635, 2015. [18] V. Noviantri, T. Tomy, and A. Chowanda. Linear and nonlinear model of brain tumor growth simulation using finite difference method. Procedia Computer Science, 179:297–304, 2021. [19] G. Powathil, M. Kohandel, S. Sivaloganathan, A. Oza, and M. Milosevic. Mathemat- ical modeling of brain tumors: effects of radiotherapy and chemotherapy. Physics in Medicine & Biology, 52(11):3291, 2007. [20] R. Rockne, EC Alvord Jr, J. K. Rockhill, and K. R. Swanson. A mathematical model for brain tumor response to radiation therapy. Journal of mathematical biology, 58:561–578, 2009. [21] R. Rockne, J. K. Rockhill, M. Mrugala, A. M. Spence, I. Kalet, K. Hendrickson, A. Lai, T. Cloughesy, EC Alvord Jr, and K. R. Swanson. Predicting the efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach. Physics in Medicine & Biology, 55(12):3271–3285, 2010. [22] B. Rysbaiuly and N. E. Mukhametkaliyeva. Iterative method of finding all thermo- physical parameters of a two-layer soil. In Journal of Physics: Conference Series, volume 2224, page 012041, 2022. [23] G. S. Stamatakos and S. G. Giatili. A numerical handling of the boundary conditions imposed by the skull on an inhomogeneous diffusion-reaction model of glioblastoma invasion into the brain: Clinical validation aspects. Cancer informatics, pages 1–16, 2017. [24] K. R. Swanson, EC Alvord Jr, and J. D. Murray. A quantitative model for differential motility of gliomas in grey and white matter. Cell proliferation, 33(5):317–329, 2000. REFERENCES 2675 [25] K. R. Swanson, EC Alvord Jr, and J. D. Murray. Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy. British journal of cancer, 86(1):14–18, 2002. [26] K. R. Swanson, C. Bridge, J. D. Murray, and EC Alvord Jr. Virtual and real brain tu- mors: using mathematical modeling to quantify glioma growth and invasion. Journal of the neurological sciences, 216(1):1–10, 2003. [27] P. Tracqui, G. C. Cruywagen, D. E. Woodward, G. T. Bartoo, J. D. Murray, and EC Alvord. A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth. Cell proliferation, 28(1):17–31, 1995. [28] M. Yamamoto and J. Zou. Simultaneous reconstruction of the initial temperature and heat radiative coefficient. Inverse problems, 17(4):1181, 2001.