# Optimal Control of Medical Drug in a Nonlocal Model of Solid Tumor Growth

Abderrahman Bouhamidi, El Harraki Imad, Yassine Melouani

► To cite this version:

Abderrahman Bouhamidi, El Harraki Imad, Yassine Melouani. Optimal Control of Medical Drug in a Nonlocal Model of Solid Tumor Growth. Journal of Mathematical Analysis and Applications, In press, 10.48550/arXiv.2503.09208. [hal-04994361](https://hal.science/hal-04994361v1)

HAL Id: hal-04994361

https://hal.science/hal-04994361v1

Submitted on 17 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal Control of Medical Drug in a Nonlocal Model of Solid Tumor Growth

Bouhamidi Abderrahman<sup>a</sup>, El Harraki Imad<sup>b</sup>, Melouani Yassine<sup>a,b</sup>

<sup>a</sup>Université du littoral côte d'opale, 50 rue Ferdinand Buisson, CS 80699, Calais, 62100, France

<sup>b</sup>Ecole nationale supérieure des mines de Rabat, Av. Hadj Ahmed Cherkaoui - B.P. : 753, Rabat, 10000, Morocco

## Abstract

This paper presents a mathematical framework for optimizing drug delivery in cancer treatment using a nonlocal model of solid tumor growth. We present a coupled system of partial differential equations that incorporate long-range cellular interactions through integral terms and drug-induced cell death. The model accounts for spatial heterogeneity in both tumor cell density and drug concentration while capturing the complex dynamics of drug resistance development. We first establish the well-posedness of the coupled system by proving the existence and uniqueness of a solution under appropriate regularity conditions. The optimal control problem is then formulated to minimize tumor size while accounting for drug toxicity constraints. Using variational methods, we derive the necessary optimality conditions and characterize the optimal control through an ad joint system. Theoretical results can help to design effective chemotherapy schedules that balance treatment efficacy with adverse effects.

**Keywords:** Optimal control, Nonlocal Tumor Growth Model, Partial Differential Equations

## Introduction

Cancer remains one of the leading causes of death worldwide, which requires continued research on tumor growth dynamics and treatment strategies. Mathematical modeling has become an essential tool for understanding complex biological processes involved and optimizing therapeutic interventions [9, 21]. The evolution of mathematical approaches to tumor growth reflects the increasing complexity of our understanding. Early models focused primarily on describing volume dynamics, with the pioneering work employing the Gompertz model [9] and von Bertalanffy model [21] establishing fundamental frameworks for characterizing growth kinetics. Although these approaches provided valuable information, the growing recognition of tumor complexity led to more sophisticated frameworks. In particular, the development of spatial models that incorporate the heterogeneous nature of the tumor microenvironment marked a significant advance. The work of Byrne and Chaplain [5] demonstrated the importance of considering spatial aspects, particularly angiogenesis, in tumor development. Building on these spatial frameworks, successive modelling efforts have included genetic and epigenetic factors to better reflect current biological understanding. Altrock et al. [1] developed sophisticated approaches that account for clonal evolution within tumors, providing critical insights into heterogeneity and treatment resistance mechanisms. This understanding was further refined by Leschiera et al. [12], whose work examining the relationship between cancer cell heterogeneity and immune response highlighted the importance of considering multiple biological scales. Although these advances significantly improved tumor growth modeling, considering only local cellular interactions overlooks interactions across various scales and distances
within the tumor environment. Such limitations may result in incomplete models that do not capture the full complexity of tumor behavior, which is crucial for understanding growth patterns and effectively optimizing treatment strategies [7]. This motivated the development of nonlocal modeling approaches capable of representing long-range cellular interactions. The groundbreaking work of Painter et al. [16] demonstrated that nonlocal terms can effectively describe pattern formation in cancer growth. Building on this foundation, Buttenschön et al. [4] established a framework for nonlocal advection models in biological systems. The contributions of Armstrong et al. [2] and Gerisch and Chaplain [8] further demonstrated the power of integro-differential equations in capturing complex tumor growth dynamics. Parallel to these developments in tumor growth modeling, the challenge of optimizing treatment protocols has led to an increasing application of control theory, particularly for chemotherapy administration. Beginning with Swan's [18] demonstration of the applicability of control theory to chemotherapy scheduling. Martin [15] extended these methods to address the crucial challenge of drug resistance, while De Pillis et al. [6] showed their utility in immunotherapy contexts. More recently, Li and You [13] extended this approach by developing optimal control strategies for reaction-diffusion models that explicitly incorporate tumor-immune interactions. The applicability of these optimal control approaches extends beyond cancer treatment to other biological systems, as demonstrated by Wu et al. [22] in their analysis of hierarchical age-structured population models in polluted environments.

In this paper, we establish a mathematical framework for analyzing nonlocal tumor growth models with chemotherapy, establishing the existence and uniqueness of a solution. Then, we develop an optimal control theory for this system, deriving necessary optimality conditions through an adjoint formulation that accounts for the nonlocal effects. Furthermore, we obtain a characterization of the optimal control strategy that provides a guidance for drug dosing protocols.

The remainder of this paper is organized as follows. In Section 1, we present the nonlocal coupled system and the optimal control framework. In Section 2, the notation and assumptions needed for our analysis are outlined. Section 3 establishes the well-posedness of the coupled system and proves the existence of an optimal control. Section 4 is devoted to deriving necessary optimality conditions and introducing an adjoint system to characterize the optimal control. In addition to the theoretical analysis, we present numerical experiments in Section 5 to validate our theoretical findings and show the efficacy of the optimal control strategy.

# 1. The Nonlocal Mathematical Model

## 1.1. Biological motivation

Cancer treatment outcomes are significantly influenced by the complex interactions between tumor cells and their microenvironment. These interactions occur not only through direct cell-to-cell contact but also via long-range signaling mechanisms that extend beyond immediate neighboring cells. Tumor cells interact through multiple mechanisms, including direct adhesion through membrane-bound molecules such as cadherins, paracrine signaling via diffusible factors, and mechanical forces transmitted through the extracellular matrix (ECM). As demonstrated by Armstrong et al. [2] and Gerisch and Chaplain [8], these interactions may be attractive or repulsive, significantly influencing both local tumor structure and overall growth patterns.

The response to treatment represents another critical aspect of tumor dynamics. Blood-borne chemotherapy agents must navigate several barriers to reach tumor cells [19]. These include diffusion through tissue, binding to target cells, and variable drug sensitivity across the tumor population. The effectiveness of treatment depends critically on both drug distribution and cellular response mechanisms, which can vary significantly across the tumor mass. Spatial heterogeneity
plays a crucial role in treatment outcomes. Tumors exhibit significant variations in cell density, drug concentration, and microenvironmental conditions throughout their structure. This heterogeneity, coupled with the complex network of cellular interactions, can lead to varying treatment responses across different regions of the tumor.

Our mathematical framework captures these biological complexities through a nonlocal velocity field $V[w_p]$. Following the approaches developed by Painter et al. [16] and Buttenschön et al. [4], we assume that cells can detect and respond to their environment within a finite range. The velocity field $V[w_p]$ integrates the effects of adhesive forces between cells, chemical gradients, and mechanical stresses over a spatial domain determined by a finite sensing radius. The model incorporates drug response through a coupled reaction-diffusion system based on the work of Trachette and Byrne [19]. This system accounts for the spatial distribution of drug concentration, cell-dependent drug uptake, and the development of resistance.

## 1.2. Model Formulation

In this subsection, we will propose a mathematical model that describes the spatiotemporal dynamics of tumor growth under the influence of chemotherapy, incorporating nonlocal effects of cell-cell interactions. Let us denote by $I$ the drug concentration in the tumor vasculature depending on the time variable $t$. The value $I(t)$ is a control variable, constrained to lie within a biologically and clinically relevant set of admissible controls. We usually, consider the set of admissible controls as:

$$ U_{ad} = \{u \in L^{\infty}(0, T) : 0 \le u(t) \le M_{tol} \text{ a.e. } t \in [0, T]\}, \quad (1) $$

where $M_{tol}$ represents the maximum tolerable drug concentration, and $T$ is any positive value. This constraint set reflects physiological limitations and safety considerations in drug administration.

The model consists of a coupled system of partial differential equations:

$$
\left\{
\begin{array}{ll}
I \in U_{ad}, & (a) \\
\displaystyle \frac{\partial p}{\partial t}(t, x) + \operatorname{div}(V[w_p(t, x)] p(t, x)) = F(p(t, x)) - C(d(t, x), p(t, x)), & t \in [0, T], x \in \mathbb{R}^N, (b) \\
\displaystyle \frac{\partial d}{\partial t}(t, x) - \operatorname{div}(D(x)\nabla d(t, x)) = \Gamma(I(t), d(t, x)), & t \in [0, T], x \in \mathbb{R}^N, (c) \\
p(0, x) = p_0(x), \quad d(0, x) = d_0(x), & x \in \mathbb{R}^N. (d) \\
\end{array}
\right.
\tag{2}
$$

Here, $p(t, x)$ represents the tumor cell density at position $x$ and time $t$, while $d(t, x)$ denotes the drug concentration. The function $V[w_p]$ is a nonlocal velocity field, where $w_p$ is a convolution functional defined as follows:

$$ w_p(t, x) = \int_{\Omega} K(x, y)p(t, y)dy, \quad (3) $$

where $K$ is a convolution kernel capturing long-range cell interactions, and $\Omega \subseteq \mathbb{R}^N$ is the domain of nonlocal sensibility, where $N \ge 1$ the spatial dimension, which is practically $N = 1, 2$. The mathematical model incorporates key biological and pharmacological mechanisms through carefully selected terms that capture the essential dynamics of tumor growth and drug treatment. The cell population dynamics are governed by two primary processes: cell proliferation, represented by the term $F(p)$, which typically follows either logistic or Gompertz growth patterns to
account for resource limitations and competition among cancer cells, and drug-induced cell death, modeled by C(d, p), which characterizes the interaction between local drug concentration and tumor cell density. The transport and metabolism of the therapeutic agent are described through two main mechanisms: First, the spatially-dependent diffusion coefficient D(x) characterizes the heterogeneous drug transport through tumor tissue, reflecting variations in tissue density, extracellular matrix composition, and interstitial pressure. Second, a nonlinear function Γ(I, d) describes the drug dynamics, incorporating various physiological processes including blood-tissue exchange, metabolic degradation, and clearance mechanisms.

From a biological and physical perspective, both the tumor cell density p(t, x) and drug concentration d(t, x) must remain non-negative throughout the temporal and spatial domains. This fundamental constraint is not automatically guaranteed by the mathematical structure of system (2). This non-negativity preservation is crucial not only for the model's physical interpretation but also for the mathematical analysis of the optimal control problem, particularly when establishing the existence of optimal solutions. In Section 3.2, we establish the non-negativity of the solutions p and d for all time t ∈ [0, T], starting with non-negative initial conditions.

### 1.3. Optimal Control Framework

To optimize the drug delivery protocol while accounting for both treatment efficacy and toxicity constraints, we formulate an optimal control problem. The objective is to minimize a cost functional that balances tumor reduction with drug-related adverse effects. Specifically, we consider the cost functional:

$$J(I) = \int_{0}^{T} \left[ \int_{\mathbb{R}^{N}} \alpha p(t, x) dx + \beta I(t)^{2} \right] dt + \gamma \int_{\mathbb{R}^{N}} p(T, x) dx. \quad (4)$$

This functional combines three key clinical considerations over the treatment period [0, T]. The first term measures the cumulative tumor density through the L¹ norm, weighted by a parameter α > 0. The second term, weighted by β > 0, penalizes high drug concentrations, reflecting the need to minimize systemic toxicity. The terminal term, weighted by γ > 0, accounts for the final tumor size, which is important for the treatment effectiveness at the conclusion of therapy [18]. Following the approaches similar to those developed by Ledzewicz and Schättler [17], the optimal control problem thus becomes the minimization of (4) subject to the control constraint system (2). Namely:

$$\left\{ \begin{array}{l} \min_{I \in U_{ad}} J(I), \\ \text{subject to} \quad (2). \end{array} \right. \qquad (5)$$

## 2. Notations and Assumptions

In this section, we will first present some notations and assumptions which will be used throughout this paper. Let define the following spaces:

Let L∞($\mathbb{R}^N$) be the space of essentially bounded measurable functions on $\mathbb{R}^N$, equipped with the norm:

$$\|f\|_{L^{\infty}(\mathbb{R}^N)} = \ess \sup_{x \in \mathbb{R}^N} |f(x)|.$$

Let C_b($\mathbb{R}^N$) be the space of continuous bounded functions on $\mathbb{R}^N$, with the uniform norm:

$$\|f\|_{C_b(\mathbb{R}^N)} = \sup_{x \in \mathbb{R}^N} |f(x)|.$$
Let $C_b^1(\mathbb{R}^N)$ be the space of continuously differentiable bounded functions on $\mathbb{R}^N$ with bounded derivatives, normed by:

$$ \|f\|_{C_b^1(\mathbb{R}^N)} = \|f\|_{C(\mathbb{R}^N)} + \|\nabla f\|_{C(\mathbb{R}^N)}. $$

For a Banach space $(X, \|\cdot\|_X)$, let $C(0, T; X)$ denote the space of continuous functions from $[0, T]$ to $X$, with the norm:

$$ \|f\|_{C(0,T;X)} = \sup_{t \in [0,T]} \|f(t)\|_X. $$

Let $L^1(0, T; X)$ be the space of Bochner integrable functions from $[0, T]$ to $X$, with the norm:

$$ \|f\|_{L^1(0,T;X)} = \int_0^T \|f(t)\|_X dt. $$

Let $W^{1,1}(\mathbb{R}^N)$ be the Sobolev space defined as:

$$ W^{1,1}(\mathbb{R}^N) = \{u \in L^1(\mathbb{R}^N) : D^\alpha u \in L^1(\mathbb{R}^N) \text{ for all } |\alpha| \le 1\}. $$

Furthermore, we state the following assumptions:

**Assumptions.** We assume the following:

(A1) $D \in C_b^1(\mathbb{R}^N)$, such that $\forall x \in \mathbb{R}^N$, $D_0 \le D(x) \le D_1$, with $D_0, D_1 > 0$,

(A2) $d_0 \in H^2(\mathbb{R}^N) \cap L^\infty(\mathbb{R}^N)$,

(A3) $\Gamma: \mathbb{R} \times \mathbb{R} \to \mathbb{R}$ is Lipschitz function w.r.t second variable,

(A4) There exists $g \in L^2(\mathbb{R}^N) \cap L^\infty(\mathbb{R}^N)$ such that:
$$ \forall (t, x) \in [0, T] \times \mathbb{R}^N, |\Gamma(I(t)), d(t, x))| \le g(x)(1 + |I(t)| + |d(t, x)|), $$

(A5) $p_0 \in L^\infty(\mathbb{R}^N) \cap W^{1,1}(\mathbb{R}^N)$,

(A6) $F: \mathbb{R} \to \mathbb{R}$ and $C: \mathbb{R} \times \mathbb{R} \to \mathbb{R}$ are Lipschitz functions (w.r.t second variable for the function $C$),

(A7) There exists $h \in L^1(\mathbb{R}^N) \cap L^\infty(\mathbb{R}^N)$ such that:
$$ \forall (t, x) \in [0, T] \times \mathbb{R}^N, |F(p(t, x)) - C(d(t, x), p(t, x))| \le h(x)(1 + |p(t, x)|), $$

(A8) $K \in C_b^2(\mathbb{R}^N \times \mathbb{R}^N) \cap L^1(\mathbb{R}^N \times \mathbb{R}^N)$, (given in (3)),

(A9) $V \in C_b^2(\mathbb{R}; \mathbb{R}^N)$.

(A10) $F \in C_b^1(\mathbb{R})$, $C \in C_b^1(\mathbb{R} \times \mathbb{R})$ and $\Gamma \in C_b^1(\mathbb{R} \times \mathbb{R})$.

(A11) $K(x, y) = 0$ for all $x \in \partial\Omega$ and $y \in \Omega$,

(A12) $V(0) = 0$ and $V'(0) = 0$.

**Remark.** In fact, the assumption (A10) will be used only in Section 4 in order to ensure the derivability of the right hand side of the system (6).
### 3. Existence and uniqueness

In this section, we establish the well-posedness of the coupled system (2). Our strategy consists of solving sequentially the parabolic equation for $d$ and then, with $d$ fixed, applying the theory developed in [10] to solve the nonlocal equation for $p$. After establishing the well-posedness of our nonlocal system, we prove the existence of an optimal control for our problem.

#### 3.1. Well-posedness of a week solution

Let us consider the coupled system (b) – (c) – (d) in (2) as:

$$ egin{cases} 
\displaystyle \frac{\partial p}{\partial t}(t, x) + \operatorname{div}(V[w_p(t, x)]p(t, x)) = F(p(t, x)) - C(d(t, x), p(t, x)), & t \in [0, T], x \in \mathbb{R}^N, & (b) \\[1.5ex] 
\displaystyle \frac{\partial d}{\partial t}(t, x) - \operatorname{div}(D(x)\nabla d(t, x)) = \Gamma(I(t), d(t, x)), & t \in [0, T], x \in \mathbb{R}^N, & (c) \\[1.5ex] 
\displaystyle p(0, x) = p_0(x), \quad d(0, x) = d_0(x), & x \in \mathbb{R}^N. & (d) 
\end{cases} 
\tag{6} $$

**Theorem 3.1** (Well-posedness of the coupled system). Let the assumptions (A1)-(A9) hold. Then for $I \in L^\infty(\mathbb{R}^N)$ fixed, the coupled system (6) has a unique solution $(p, d)$ in $C(0, T; L^1(\mathbb{R}^N)) \times C(0, T; H^1(\mathbb{R}^N))$. Furthermore, we have the following estimates:

$$ \|p(t, \cdot)\|_{L^1(\mathbb{R}^N)} \le C_1 (\|p_0\|_{L^1(\mathbb{R}^N)} + T), \tag{7} $$

$$ \|p(t, \cdot)\|_{L^{\infty}(\mathbb{R}^N)} \le C_2 (\|p_0\|_{L^{\infty}(\mathbb{R}^N)} + T), \tag{8} $$

$$ \|d(t, \cdot)\|_{L^2(\mathbb{R}^N)} \le C_3 (1 + \|d_0\|_{L^2(\mathbb{R}^N)} + \|I\|_{L^2(0,T)}), \tag{9} $$

and

$$ \|d\|_{L^{\infty}(0,T;L^{\infty}(\mathbb{R}^N))} \le C_4 (1 + \|d_0\|_{L^{\infty}(\mathbb{R}^N)} + \|I\|_{L^{\infty}(0,T)}). \tag{10} $$

where $C_1, C_2, C_3$, and $C_4$ are positive constants.

*Proof.* • We start by proving the well-posedness of the parabolic equation (6)-(c).

Let $\mathcal{A}d = -\operatorname{div}(D(x)\nabla d)$ with domain $D(\mathcal{A}) = H^2(\mathbb{R}^N)$. Under Assumption (A1), this operator is uniformly elliptic. By Theorem 2.5.1 in Zheng [23], under Assumptions (A1)-(A3), there exists a unique solution $d \in C(0, T; H^1(\mathbb{R}^N))$ which can be represented as follows:

$$ d(t, x) = e^{-t\mathcal{A}}d_0 + \int_0^t e^{-(t-s)\mathcal{A}}\Gamma(I(s), d(s, \cdot))ds. \tag{11} $$

For the $L^2$ estimate, we first note that $\{e^{-t\mathcal{A}}\}_{t \ge 0}$ is a strongly continuous semigroup on $L^2(\mathbb{R}^N)$ with the property:

$$ \|e^{-t\mathcal{A}}\|_{L^2 \to L^2} \le 1, $$
due to the uniform ellipticity of $\mathcal{A}$. By taking the $L^2$ norm of the mild solution we get:

$$ \begin{aligned} \|d(t, \cdot)\|_{L^2(\mathbb{R}^N)} &\le \|e^{-t\mathcal{A}}d_0\|_{L^2(\mathbb{R}^N)} + \int_0^t \|e^{-(t-s)\mathcal{A}}\Gamma(I(s), d(s, \cdot))\|_{L^2(\mathbb{R}^N)}ds \\ &\le \|d_0\|_{L^2(\mathbb{R}^N)} + \int_0^t \|\Gamma(I(s), d(s, \cdot))\|_{L^2(\mathbb{R}^N)}ds \end{aligned} $$

By Assumption (A4), we have:

$$ \|\Gamma(I(s), d(s, \cdot))\|_{L^2(\mathbb{R}^N)} \le M_1(1 + \|I\|_{L^\infty(0,T)} + \|d(s, \cdot)\|_{L^2(\mathbb{R}^N)}) $$

By inserting this into our previous inequality we have:

$$ \begin{aligned} \|d(t, \cdot)\|_{L^2(\mathbb{R}^N)} &\le \|d_0\|_{L^2(\mathbb{R}^N)} + M_1 \int_0^t (1 + \|I\|_{L^\infty(0,T)} + \|d(s, \cdot)\|_{L^2(\mathbb{R}^N)})ds \\ &\le \|d_0\|_{L^2(\mathbb{R}^N)} + M_1 t(1 + \|I\|_{L^\infty(0,T)}) + M_1 \int_0^t \|d(s, \cdot)\|_{L^2(\mathbb{R}^N)}ds \end{aligned} $$

Therefore, applying Gronwall's inequality leads to:

$$ \|d(t, \cdot)\|_{L^2(\mathbb{R}^N)} \le (\|d_0\|_{L^2(\mathbb{R}^N)} + M_1 T(1 + \|I\|_{L^\infty(0,T)}))e^{M_1 T} $$

This gives the estimate (9). For the $L^\infty$ estimate, we use that the heat semigroup satisfies the maximum principle:

$$ \|e^{-t\mathcal{A}}\|_{L^\infty \to L^\infty} \le 1, $$

Then:

$$ \begin{aligned} \|d(t, \cdot)\|_{L^\infty(\mathbb{R}^N)} &\le \|d_0\|_{L^\infty(\mathbb{R}^N)} + \int_0^t \|\Gamma(I(s), d(s, \cdot))\|_{L^\infty(\mathbb{R}^N)}ds \\ &\le \|d_0\|_{L^\infty(\mathbb{R}^N)} + M_2 \int_0^t (1 + \|I\|_{L^\infty(0,T)} + \|d(s, \cdot)\|_{L^\infty(\mathbb{R}^N)})ds \end{aligned} $$

Therefore, applying Gronwall's inequality yields the estimate (10).

* Now, for the well-posedness of the nonlocal equation (6)-(b), let $d \in C(0, T; H^1(\mathbb{R}^N)) \cap L^\infty((0, T); L^\infty(\mathbb{R}^N))$ be fixed, then we will use Banach's fixed theorem and Keimer et al.[11] to establish the desired result. In fact, Let $M := 2\|p_0\|_{L^1(\mathbb{R}^N)}$, we define the following closed ball subset:

$$ E = \{p \in C([0, T]; L^1(\mathbb{R}^N)) : \|p(t, \cdot)\|_{L^1(\mathbb{R}^N)} \le M\}, \quad (12) $$

of the Banach space $C([0, T]; L^1(\mathbb{R}^N))$. Furthermore, we consider the map $\Phi: E \to C([0, T]; L^1(\mathbb{R}^N))$ defined by:

$$ \Phi(p)(t, x) = p_0(X[t, x](0)) \det(D_x X[t, x](0)) + \int_0^t [F(p) - C(d, p)](s, X[t, x](s)) \det(DX[t, x](s))ds, \quad (13) $$
where $X[t, x]$ solves the characteristic equation:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial X[t, x](s)}{\partial s} = V[w_p(s, X[t, x](s))], \quad s \in [0, T], \\ X[t, x](t) = x, \end{array} \right. \qquad (14) $$

and $V$ is given in (2). By using the characteristics, for $p \in E$, and for any $t \in [0, T]$ we obtain:

$$ \begin{aligned} \|\Phi(p)(t, \cdot)\|_{L^1(\mathbb{R}^N)} &\le \int_{\mathbb{R}^N} |p_0(X[t, x](0))| \det(D_x X[t, x](0)) dx \\ &\quad + \int_{\mathbb{R}^N} \int_0^t |[F(p) - C(d, p)](\tau, X[t, x](\tau))| \det(D_x X[t, x](\tau)) d\tau dx. \end{aligned} $$

Now, the change of variables $y = X[t, x](0)$ for the first term and $y = X[t, x](\tau)$ for the second leads to:

$$ \int_{\mathbb{R}^N} |p_0(X[t, x](0))| \det(D_x X[t, x](0)) dx = \|p_0\|_{L^1(\mathbb{R}^N)}. $$

For the second term by using the assumption (A7), we obtain:

$$ \begin{aligned} \|\Phi(p)(t, \cdot)\|_{L^1(\mathbb{R}^N)} &\le \|p_0\|_{L^1(\mathbb{R}^N)} + T \int_{\mathbb{R}^N} h(x)(1 + |p(t, x)|)dx \\ &\le \|p_0\|_{L^1(\mathbb{R}^N)} + C_1 T(1 + M), \end{aligned} $$

where $C_1 = \max \{\|h\|_{L^1(\mathbb{R}^N)}, \|h\|_{L^\infty(\mathbb{R}^N)}\}$. By choosing $T$ small enough so that $C_1T(1+M) \le \|p_0\|_{L^1(\mathbb{R}^N)}$, we get for $t \in [0, T]$ the following:

$$ \|\Phi(p)(t, \cdot)\|_{L^1(\mathbb{R}^N)} \le 2\|p_0\|_{L^1(\mathbb{R}^N)} = M. $$

Then $\Phi(E) \subseteq E$. As $F$ and $C$ are Lipschitz functions with the Lipschitz constants $L_F$ and $L_C$ respectively. Then for $p_1, p_2 \in E$ and for $t \in [0, T]$, we have:

$$ \begin{aligned} \|\Phi(p_1)(t, \cdot) - \Phi(p_2)(t, \cdot)\|_{L^1(\mathbb{R}^N)} &= \int_{\mathbb{R}^N} \int_0^t |F(p_1) - F(p_2) - (C(d, p_1) - C(d, p_2))| \cdot \det(D_x X[t, x](\tau)) d\tau dx \\ &\le (L_F + L_C) \int_0^t \|p_1(\tau, \cdot) - p_2(\tau, \cdot)\|_{L^1(\mathbb{R}^N)} d\tau. \end{aligned} $$

Therefore:

$$ \|\Phi(p_1) - \Phi(p_2)\|_{C([0,T];L^1(\mathbb{R}^N))} \le (L_F + L_C)T \|p_1 - p_2\|_{C([0,T];L^1(\mathbb{R}^N))}. $$

Then, by choosing $T^* = \min \left\{ T, \frac{1}{2(L_F + L_C)} \right\}$, $\Phi$ is a contraction on $E$ for all $t$ in $[0, T^*]$. By using Banach's fixed point theorem, there exists a unique $p^* \in E$ such that $\Phi(p^*) = p^*$.

Consequently, using assumptions (A1)-(A9) and according to Theorem 3.24 in Keimer et al. [11], we establish a unique solution for the nonlocal equation (6) in $C([0, T^*]; L^1(\mathbb{R}^N))$, furthermore, the existence and uniqueness of the solution holds for any $T \ge 0$.
To obtain the estimates for $p$, we use the characteristics representation as follows, for all $t \ge t_0$ where $t_0 \in [0, T]$:

$$p(t, x) = p(t_0, X[t, x](t_0)) \det(D_x X[t, x](t_0)) + \int_{t_0}^{t} [F(p) - C(d, p)](\tau, X[t, x](\tau)) \det(D_x X[t, x](\tau)) d\tau. \quad (15)$$

By taking the $L^1$ norm and using the change of variables $y = X[t, x](0)$ for the first term and $y = X[t, x](\tau)$ for the second we get the following:

$$\|p(t, \cdot)\|_{L^1(\mathbb{R}^N)} \le \|p_0\|_{L^1(\mathbb{R}^N)} + \int_0^t \|F(p(\tau)) - C(d(\tau), p(\tau))\|_{L^1(\mathbb{R}^N)} d\tau.$$

By using (A7), we get:

$$\|p(t, \cdot)\|_{L^1(\mathbb{R}^N)} \le \|p_0\|_{L^1(\mathbb{R}^N)} + C_1 \int_0^t (1 + \|p(\tau)\|_{L^1(\mathbb{R}^N)}) d\tau.$$

Therefore, by using Gronwall's inequality, we get the desired $L^1$ estimate (7). Similarly, to get the $L^\infty$ estimate we have:

$$|p(t, x)| \le |p_0(\xi[t, x](0))| \cdot \det(D_x X[t, x](0)) + \int_0^t |F(p) - C(d, p)| (\tau, X[t, x](\tau)) \cdot \det(D_x X[t, x](\tau)) d\tau.$$

From Lemma 3.3 in [11], we have:

$$\det(D_x X[t, x](\tau)) = \exp \left( \int_{\tau}^{t} - \operatorname{div} \left( V(w(s, X[t, x](s))) \right) ds \right). \quad (16)$$

According to assumption (A9), we have:

$$|\det(D_x X[t, x](\tau))| \le \exp(LT),$$

where $L$ depends on bounds of derivatives of $V$. Therefore:

$$|p(t, x)| \le \|p_0\|_{L^\infty(\mathbb{R}^N)} \exp(LT) + \int_0^t \exp(LT) |F(p) - C(d, p)| (\tau, X[t, x](\tau)) d\tau,$$

by taking the supremum over $x \in \mathbb{R}^N$, and using Gronwall's inequality, we get the desired $L^\infty$ estimate (8). $\square$

## 3.2. Existence of an Optimal Control

Having established the well-posedness of our coupled system, we now turn our attention to the existence of an optimal control. We consider the following optimal control problem consisting on the minimization of the cost functional (4) subject to the system (2).

Before announcing our main result in this section, we need to establish the non-negativity of $p$ as given in the following lemma:

**Lemma 3.1 (Non-negativity of the solutions).** For $I \in \mathcal{U}_{ad}$ fixed, let the assumptions (A1)-(A9) hold. If the initial conditions are non-negative $d_0 \ge 0$ and $p_0 \ge 0$, then for any $(t, x) \in [0, T] \times \mathbb{R}^N$, $p(t, x)$ and $d(t, x)$ are non-negative.
Proof. Let us first establish the non-negativity of $d$. Let $\mathcal{A}d = -\operatorname{div}(D(x)\nabla d)$. By the parabolic maximum principle, the semigroup $e^{-t\mathcal{A}}_{t\ge 0}$ is positivity-preserving: if $f \ge 0$, then $e^{-t\mathcal{A}}f \ge 0$ for all $t \ge 0$. We will proceed by contradiction. Let $t^*$ be the supremum of all times $t \in [0, T]$ such that $d(s, x) \ge 0$ for all $s \in [0, t]$ and $x \in \mathbb{R}^N$. Since $d_0 \ge 0$ and $d$ is continuous, $t^*$ is well-defined and positive. Suppose by contradiction that $t^* < T$. By continuity of $d$, we have $d(t^*, x) \ge 0$ for all $x \in \mathbb{R}^N$, and by maximality of $t^*$, there exist $\varepsilon > 0$ and $x^* \in \mathbb{R}^N$ such that $d(t, x^*) < 0$ for some $t \in [t^*, t^* + \varepsilon]$. For $t \in [t^*, t^* + \varepsilon]$, we use the variation of constants formula:

$$d(t, x) = e^{-(t-t^*)\mathcal{A}}d(t^*, \cdot) + \int_{t^*}^{t} e^{-(t-s)\mathcal{A}}\Gamma(I(s), d(s, \cdot))ds. \quad (17)$$

Let $d^-(t, \cdot) = \max(-d(t, \cdot), 0)$ denote the negative part of $d$. Since $\Gamma$ is Lipschitz in the second variable with constant $L$, we have:

$$\Gamma(I(s), d(s, \cdot)) = \Gamma(I(s), 0) + \int_{0}^{d(s, \cdot)} \phi(I(s), \tau)d\tau$$

where $|\phi(I(s), \tau)| \le L$ is the derivative of $\Gamma$ with respect to the second variable (which exists almost everywhere by Rademacher's theorem since $\Gamma$ is Lipschitz). By the growth condition we have $|\Gamma(I(s), 0)| \le g(x)(1 + |I(s)|)$, therefore:

$$\Gamma(I(s), d(s, x)) \ge -g(x)(1 + |I(s)|) - L|d(s, x)|.$$

By using the positivity-preserving property of the semigroup and the non-negativity of $d(t^*, \cdot)$ we get:

$$\begin{aligned} d(t, x) &\ge -\int_{t^*}^{t} e^{-(t-s)\mathcal{A}}[g(\cdot)(1 + |I(s)|) + L|d(s, \cdot)|]ds \\ &\ge -C_1 \int_{t^*}^{t} (1 + \|d(s, \cdot)\|_{L^\infty(\mathbb{R}^N)})ds, \end{aligned}$$

where $C$ depends on $\|g\|_{L^\infty(\mathbb{R}^N)}$, $\|I\|_{L^\infty(0,T)}$, and $L$. Therefore:

$$\|d^-(t, \cdot)\|_{L^\infty(\mathbb{R}^N)} \le C_1 \int_{t^*}^{t} (1 + \|d(s, \cdot)\|_{L^\infty(\mathbb{R}^N)})ds. \quad (18)$$

Now, in the region where $d$ is negative (i.e., for $t \in [t^*, t^* + \varepsilon]$), we have $d(t, x) = -d^-(t, x)$, and therefore:

$$\|d(s, \cdot)\|_{L^\infty(\mathbb{R}^N)} = \|d^-(s, \cdot)\|_{L^\infty(\mathbb{R}^N)}.$$

By Substituting this back into (18) and applying Gronwall's inequality we get the following:

$$\|d^-(t, \cdot)\|_{L^\infty(\mathbb{R}^N)} \le C_1(t - t^*) \exp(C_1(t - t^*)). \quad (19)$$

For sufficiently small $\varepsilon > 0$, this implies $d^-(t, \cdot) = 0$ for $t \in [t^*, t^* + \varepsilon]$, contradicting our assumption that $d$ becomes strictly negative after $t^*$. Therefore, $t^* = T$ and $d$ remains non-negative on $[0, T]$.

Now we will prove the non-negativity of $p$. For $t \in [0, T]$, let $p^+(t, \cdot)$ and $p^-(t, \cdot)$ denote the usual positive and negative parts of $p(t, \cdot)$, respectively. We have: $p^+(t, \cdot) := \max(p(t, \cdot), 0) \ge 0$,
and $p^-(t, \cdot) := \max(-p(t, \cdot), 0) = -\min(p(t, \cdot), 0) \ge 0$ and $p(t, \cdot) = p^+(t, \cdot) - p^-(t, \cdot)$. For any fixed $t^* \in [0, T]$ we have:

$$ \|p^-(t^*, \cdot)\|_{L^1(\mathbb{R}^N)} = \int_{\mathbb{R}^N} |p^-(t^*, x)|dx = - \int_{\{p(t^*, \cdot) < 0\}} p(t^*, x)dx. $$

By using the characteristic formula (15) and integrating on $\{p(t^*, \cdot) < 0\}$, we get the following:

$$ \begin{aligned} \|p^-(t^*, \cdot)\|_{L^1(\mathbb{R}^N)} = & - \int_{\{p(t^*, \cdot) < 0\}} p_0(X[t^*, x](0)) \det(D_x X[t^*, x](0)) dx \\ & - \int_{\{p(t^*, \cdot) < 0\}} \int_0^{t^*} [F(p) - C(d, p)](\tau, X[t^*, x](\tau)) \det(D_x X[t^*, x](\tau)) d\tau dx. \end{aligned} $$

Since $p_0 \ge 0$ and according to (16), the determinant is non negative, then:

$$ - \int_{\{p(t^*, \cdot) < 0\}} p_0(X[t^*, x](0)) \det(D_x X[t^*, x](0)) dx \le 0. $$

By using the change of variables $y = X[t^*, x](\tau)$, the assumption (A6), and the estimate (7) respectively, we obtain:

$$ \|p^-(t^*, \cdot)\|_{L^1(\mathbb{R}^N)} \le \int_{\mathbb{R}^N} \int_0^{t^*} \left| [F(p) - C(d, p)](\tau, X[t^*, x](\tau)) \right| \det(D_x X[t^*, x](\tau)) d\tau dx. $$

Using the assumption (A6) and the estimate (7), we get:

$$ \|p^-(t^*, \cdot)\|_{L^1(\mathbb{R}^N)} \le C_1 t^* (1 + C_2 (\|p_0\|_{L^1(\mathbb{R}^N)} + T)) \le C' t^*, $$

where $C_1, C_1$, and $C'$ are positive constants. Let $\varepsilon > 0$, by choosing $t^* > 0$ small enough such that $C't^* < \varepsilon$, we get $\|p^-(t)\|_{L^1(\mathbb{R}^N)} \le C't^* < \varepsilon$, for all $t \in [0, t^*]$. It follows that $p^-(t, \cdot) = 0$ a.e on $\mathbb{R}^N$, for all $t \in [0, t^*]$. Finally, $p(t, \cdot) \ge 0$ a.e on $\mathbb{R}^N$, for all $t \in [0, t^*]$. The continuity of $p$ ensures that $p(t, \cdot) \ge 0$ on $\mathbb{R}^N$, for all $t \in [0, t^*]$.

Now, for all $t \in [t^*, 2t^*]$, by using (15) we get similarly:

$$ p(t, x) = p(t^*, X[t, x](t^*)) + \int_{t^*}^{t} [F(p) - C(d, p)](\tau, X[t, x](\tau)) \det(D_x X[t, x](\tau)) d\tau. \quad (20) $$

Since $p(t^*, x) \ge 0$ for all $x \in \mathbb{R}^N$, by following the same arguments as in the interval $[0, t^*]$, we get:

$$ \|p^-(t, \cdot)\|_{L^1(\mathbb{R}^N)} \le C'(t - t^*) \le C' t^* < \varepsilon. $$

Therefore $p \ge 0$ on $[t^*, 2t^*]$. Let $L$ be the greatest integer in $\mathbb{N}$ such that $Lt^* \le T$. We can iterate the argument on each interval $[kt^*, (k+1)t^*]$ using the fact that; $p(kt^*) \ge 0$ as initial condition at each step for $k = 0, 1, \dots, L-1$, and on $[Lt^*, T]$ using $p(Lt^*) \ge 0$ as initial condition. Therefore we obtain the desired result $p(t, \cdot) \ge 0$ on $\mathbb{R}^N$ for all $t \in [0, T]$. $\square$

Now, we can establish the result of existence of an optimal control:

**Theorem 3.2 (Existence of an Optimal Control).** Let the assumptions (A1)-(A9) hold. Then there exists an optimal control $I^* \in \mathcal{U}_{ad}$ that minimizes the cost functional $J$.
Proof. First, we establish that the cost functional $J$ is bounded below. Due to the non-negativity of $p$ (see Lemma 3.1), we have $J(I) \ge 0$ for all $I \in \mathcal{U}_{ad}$. This ensures that $\inf_{I \in \mathcal{U}_{ad}} J(I)$ is well-defined. Let $(I^n)_n \subset \mathcal{U}_{ad}$ be a minimizing sequence, i.e.,

$$ \lim_{n \to \infty} J(I^n) = \inf_{I \in \mathcal{U}_{ad}} J(I). \qquad (21) $$

For each $n$, let $(p^n, d^n)$ be the solution to system (6) corresponding to the control $I^n$. From estimates (7),(8) and (9) in Theorem 3.1, we obtain the uniform bounds:

$$ \|p^n\|_{L^\infty(0,T;L^1(\mathbb{R}^N))} \le C_1, \quad \|p^n\|_{L^\infty(0,T;L^\infty(\mathbb{R}^N))} \le C_2, \quad \|d^n\|_{L^2(0,T;L^2(\mathbb{R}^N))} \le C_3. \qquad (22) $$

where $C_1, C_2$, and $C_3$ are some positive constants, and we have also $\|I^n\|_{L^\infty(0,T)} \le M_{tol}$ where $M_{tol}$ is defined in (1). By the Banach-Alaoglu theorem [3], there exist subsequences (still denoted by the same notations) and limit functions such that:

$$ I^n \rightharpoonup I^* \quad \text{weakly}^* \text{ in } L^\infty(0,T), \qquad (23) $$

$$ p^n \rightharpoonup p^* \quad \text{weakly}^* \text{ in } L^\infty(0,T; L^\infty(\mathbb{R}^N)), \qquad (24) $$

$$ d^n \rightharpoonup d^* \quad \text{weakly} \text{ in } L^2(0,T; L^2(\mathbb{R}^N)). \qquad (25) $$

Now we will show that $(p^*, d^*)$ is the solution of the system (6)-(b)-(c)-(d) for $I = I^*$. Indeed, to show that $p^*$ and $d^*$ satisfy the initial conditions (d)-(6), we consider a test function $\phi \in C_c^1([-T, T] \times \mathbb{R}^N)$, then integrating by part with respect to $t$:

$$ \int_0^T \int_{\mathbb{R}^N} \partial_t d^n(t, x) \phi(t, x) dxdt = \int_{\mathbb{R}^N} d^n(0, x) \phi(0, x) dx - \int_0^T \int_{\mathbb{R}^N} d^n(t, x) \partial_t \phi(t, x) dxdt. $$

According to weak convergence (25) and passing to the limit we get:

$$ \int_0^T \int_{\mathbb{R}^N} \partial_t d^*(t, x) \phi(t, x) dxdt = \lim_{n \to \infty} \int_{\mathbb{R}^N} d^*(0, x) \phi(0, x) dx - \int_0^T \int_{\mathbb{R}^N} d^*(t, x) \partial_t \phi(t, x) dxdt. $$

Integrating again by part the second integral in the right hand side of last relation, we get:

$$ \lim_{n \to \infty} \int_{\mathbb{R}^N} [d^n(0, x) - d^*(0, x)] \phi(0, x) dx = 0. $$

As $d^n(0, x) = d_0(x)$, for all $n \ge 0$, then

$$ \int_{\mathbb{R}^N} [d_0(x) - d^*(0, x)] \phi(0, x) dx = 0, $$

which shows that $d^*$ satisfies the initial condition. By using the same argument we have also $p^*$ satisfying the initial condition (d)-(6). Now it remains to show that $(p^*, d^*)$ is the solution of the system (6)-(b)-(c). For all $\phi \in C_c^1([0, T] \times \mathbb{R}^N)$, from equation (c)-(6), we have:

$$ \int_0^T \int_{\mathbb{R}^N} \partial_t d^n \phi \, dxdt + \int_0^T \int_{\mathbb{R}^N} D(x) \nabla d^n \cdot \nabla \phi \, dxdt = \int_0^T \int_{\mathbb{R}^N} \Gamma(I^n, d^n) \phi \, dxdt $$

Using the weak* convergence (23), the weak convergence (25), and the continuity of $\Gamma$ we obtain the limit equation:

$$ \int_0^T \int_{\mathbb{R}^N} \partial_t d^* \phi \, dxdt + \int_0^T \int_{\mathbb{R}^N} D(x) \nabla d^* \cdot \nabla \phi \, dxdt = \int_0^T \int_{\mathbb{R}^N} \Gamma(I^*, d^*) \phi \, dxdt. \qquad (26) $$
For the nonlocal equation (b)-(6), for all $\phi \in C_c^1([0, T] \times \mathbb{R}^N)$, we have:

$$ \int_0^T \int_{\mathbb{R}^N} \partial_t p^n \phi \, dxdt \to \int_0^T \int_{\mathbb{R}^N} \partial_t p^* \phi \, dxdt. \quad (27) $$

For the nonlocal term $w_{p^n}$, by using Fubini's theorem we have:

$$ \begin{aligned} \int_0^T \int_{\mathbb{R}^N} w_{p^n}(t, x)\phi(t, x)dxdt &= \int_0^T \int_{\mathbb{R}^N} \int_{\Omega} K(x, y)p^n(t, y)\phi(t, x)dydxdt, \\ &= \int_0^T \int_{\Omega} p^n(t, y)\left(\int_{\mathbb{R}^N} K(x, y)\phi(t, x)dx\right)dydt. \end{aligned} $$

Let define $\psi(t, y) = \int_{\mathbb{R}^N} K(x, y)\phi(t, x)dx$, by assumption (A8), we have

$$ \int_0^T \int_{\mathbb{R}^N} |\psi(t, y)|dydt < \infty. $$

Therefore $\psi \in L^1(0, T; L^1(\Omega))$, and by weak* convergence of $p^n$ in $L^\infty(0, T; L^\infty(\Omega))$ (which follows from the weak* convergence in $L^\infty(0, T; L^\infty(\mathbb{R}^N))$). It follows that:

$$ \int_0^T \int_{\Omega} p^n(t, y)\psi(t, y)dydt \to \int_0^T \int_{\Omega} p^*(t, y)\psi(t, y)dydt, $$

which shows that $w_{p^n} \rightharpoonup w_{p^*}$ weakly* in $L^\infty(0, T; L^\infty(\mathbb{R}^N))$. From assumption (A9), the continuity of V ensures that $V[w_{p^n}] \rightharpoonup V[w_{p^*}]$ weakly* in $L^\infty(0, T; L^\infty(\mathbb{R}^N))$. Now, we have:

$$ \int_0^T \int_{\mathbb{R}^N} p^n V[w_{p^n}] \cdot \nabla \phi \, dxdt - \int_0^T \int_{\mathbb{R}^N} p^* V[w_{p^*}] \cdot \nabla \phi \, dxdt \quad (28) $$

$$ = \int_0^T \int_{\mathbb{R}^N} (p^n - p^*) V[w_{p^n}] \cdot \nabla \phi \, dxdt + \int_0^T \int_{\mathbb{R}^N} p^*(V[w_{p^n}] - V[w_{p^*}]) \cdot \nabla \phi \, dxdt. \quad (29) $$

For the first term we have $V[w_{p^n}] \cdot \nabla \phi \in L^1(0, T; L^1(\mathbb{R}^N))$ since V is bounded, hence by the weakly convergence* (24) we get:

$$ \int_0^T \int_{\mathbb{R}^N} (p^n - p^*) V[w_{p^n}] \cdot \nabla \phi \, dxdt \to 0. \quad (30) $$

For the second term, since $V[w_{p^n}] \rightharpoonup V[w_{p^*}]$ weakly* in $L^\infty(0, T; L^\infty(\mathbb{R}^N))$ and $p^* \in L^\infty(0, T; L^\infty(\mathbb{R}^N))$ we have:

$$ \int_0^T \int_{\mathbb{R}^N} p^*(V[w_{p^n}] - V[w_{p^*}]) \cdot \nabla \phi \, dxdt \to 0. \quad (31) $$

Hence:

$$ \int_0^T \int_{\mathbb{R}^N} p^n V[w_{p^n}] \cdot \nabla \phi \, dxdt \to \int_0^T \int_{\mathbb{R}^N} p^* V[w_{p^*}] \cdot \nabla \phi \, dxdt. \quad (32) $$
For the right hand side, the use of the weak* convergence and the continuity of F and C (assumption (A6)) allows passage to the limit:

$$ \int_{0}^{T} \int_{\mathbb{R}^{N}} (F(p^{n}) - C(d^{n}, p^{n})) \phi \, dxdt \rightarrow \int_{0}^{T} \int_{\mathbb{R}^{N}} (F(p^{*}) - C(d^{*}, p^{*})) \phi \, dxdt. \quad (33) $$

From (26), (27), (32), and (33), we conclude that $p^*$ and $d^*$ satisfy the system (b)-(c)(6).

Now we will show that the cost functional J is lower semicontinuous. For the space-time integral term, we use the non-negativity of p (from 3.1) and monotone convergence theorem to handle the integration over $\mathbb{R}^N$. This, combined with the weak* convergence of $p^n$ on bounded domains, leads to:

$$ \int_{0}^{T} \int_{\mathbb{R}^{N}} \alpha p^{*}(t, x) dxdt \leq \liminf_{n \to \infty} \int_{0}^{T} \int_{\mathbb{R}^{N}} \alpha p^{n}(t, x) dxdt. $$

The quadratic control term is convex, hence weakly lower semicontinuous, giving:

$$ \int_{0}^{T} \beta(I^{*}(t))^{2} dt \leq \liminf_{n \to \infty} \int_{0}^{T} \beta(I^{n}(t))^{2} dt. $$

The terminal term follows similarly to the first term. Therefore by combining these inequalities we get:

$$ J(I^*) \leq \liminf_{n \to \infty} J(I^n) = \inf_{I \in \mathcal{U}_{ad}} J(I), $$

which proves the optimality of $I^*$.

# 4. Optimality Conditions

In this section, we will derive the necessary optimality conditions for our optimal control problem. Our approach follows the classical framework of optimal control theory for partial differential equations, as outlined in the seminal works of Lions [14] and Tröltzsch [20]. The derivation of these conditions typically involves three main steps, starting by proving the differentiability of the control-to-state mapping, introducing an adjoint system, and deriving the optimality conditions.

## 4.1. Differentiability of the Control-to-State Mapping

Let us define the following control-to-state mapping:

$$ \left\{ \begin{array}{ll} G: & \mathcal{U}_{ad} \longrightarrow C(0, T; L^1(\mathbb{R}^N)) \times C([0, T]; H^1(\mathbb{R}^N)), \\ I & \longrightarrow G(I) = (p, d). \end{array} \right. \qquad (34) $$

where $(p, d)$ is the solution to system (6) corresponding to a given $I \in U_{ad}$. Let us consider the notation $\frac{\partial G}{\partial h}(I)$ to denote the gâteaux derivative of $G$ at $I \in U_{ad}$ in the direction $h \in L^\infty(0, T)$. We have the following lemma:

**Lemma 4.1.** Let Assumptions (A1)-(A10) hold. Then for any $h \in L^\infty(0, T)$, the linearized system:

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial \xi}{\partial t}(t, x) + \operatorname{div}((V[w_p]\xi)(t, x)) + \operatorname{div}((V'[w_p]w_\xi p)(t, x)) = g(d, p)(t, x)\xi(t, x) - \partial_1 C(d, p)(t, x)\eta(t, x), & (a) \\ \\ \displaystyle \frac{\partial \eta}{\partial t}(t, x) - \operatorname{div}(D(x)\nabla\eta(t, x)) = \partial_1 \Gamma(I, d)(t, x)h(t) + \partial_2 \Gamma(I, d)(t, x)\eta(t, x), & (b) \\ \\ \xi(0, x) = 0, \quad \eta(0, x) = 0, \end{array} \right. $$
has a unique solution $(\xi, \eta)$ in $C(0, T; L^1(\mathbb{R}^N)) \times C(0, T; H^1(\mathbb{R}^N))$. Where $g(d,p) := F'(p) - \partial_2 C(d,p)$. $V'$, $F'$ are the first derivatives of $V$ and $F$ respectively, while $\partial_1$ and $\partial_2$ denote the first derivative with respect to first and second variable, respectively.

**Proof.** The second equation (35)-(b) is a linear parabolic equation, then it has a unique solution. For the well-posedness of equation (35)-(a), the proof is based on the work of Keimer et al. [11]. However, our equation (35)-(a) still slightly different from the one discussed in [11]. Let us consider the following equation:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \xi}{\partial t}(t, x) + \operatorname{div}((V[w_p]\xi)(t, x)) = -\operatorname{div}((V'[w_p]\operatorname{wp})(t, x)) - \partial_1 C(d,p)(t, x)\eta(t, x) - g(d,p)(t, x)\xi(t, x), \\ \xi(0, x) = 0, \end{array} \right. \qquad (36) $$

where $w$ is a fixed function in $C(0, T; C_b^1(\mathbb{R}^N))$ independent from $\xi$. The equation (36) is a linear balance hyperbolic law. Hence, by using the Lemma 3.4 in [11] and under the assumptions (A1)-(A9), we deduce that (36) has a unique solution which may be written as follows:

$$ \xi(t, x) = - \int_{0}^{t} \left[ \partial_{1} C(d, p)\eta + \operatorname{div}(V'[w_p] \operatorname{wp}) \right] (s, X[t, x](s)) \exp \left( \int_{s}^{t} g(d, p)(\tau, X[t, x](\tau)) d\tau \right) \det(D_x X[t, x](s)) ds, \qquad (37) $$

where $X[t, x]$ is the solution of the characteristic equation (14). Now let us consider the following closed ball:

$$ B = \{ f \in C(0, T; C_b^1(\mathbb{R}^N)) : \|f\|_{C(0, T; C_b^1(\mathbb{R}^N))} \le M \}, \qquad (38) $$

where $M$ is a positive constant. Let consider the map $L : B \to B$ defined by:

$$ L(w)(t, x) = - \int_{0}^{t} \int_{X[t, \Omega](\tau)} K(x, X[\tau, y](t)) \left[ \partial_{1} C(d, p)\eta + \operatorname{div}(V'[w_p] \operatorname{wp}) \right] (\tau, y) \exp \left( \int_{s}^{t} g(d, p)(\tau, X[\tau, y](s)) ds \right) dy d\tau, \qquad (39) $$

where $\Omega \subseteq \mathbb{R}^N$ is the domain appearing in (3), and $X[t, \Omega](\tau)$ is the image by $X[t, \cdot](\tau)$ of $\Omega$. We first show that the mapping $L$ is well-defined, for $w \in B$ we have:

$$ \begin{aligned} |L(w)(t, x)| & \le \int_{0}^{t} \int_{\mathbb{R}^N} \left| K(x, X[\tau, y](t)) (\partial_1 C(d, p)\eta)(\tau, y) \exp \left( \int_{s}^{t} g(d, p)(\tau, X[\tau, y](s)) ds \right) \right| dy d\tau \\ & \quad + \int_{0}^{t} \int_{\mathbb{R}^N} \left| K(x, X[\tau, y](t)) \operatorname{div}(V'[w_p] \operatorname{wp})(\tau, y) \exp \left( \int_{s}^{t} g(d, p)(\tau, X[\tau, y](s)) ds \right) \right| dy d\tau. \end{aligned} $$

According to the higher regularity Theorem 4.2 in [11] and the assumptions (A5), (A8), (A9) and (A10), we obtain that $p \in C(0, T; W^{1,1}(\mathbb{R}^N))$. Then

$$ A_1 := \|K\|_{L^\infty(\mathbb{R}^N \times \mathbb{R}^N)} \|\partial_1 C(d,p)\|_{L^\infty(0,T,L^\infty(\mathbb{R}^N))} \|\eta\|_{L^\infty(0,T,L^\infty(\mathbb{R}^N))} \exp \left( \|g(d,p)\|_{L^1(0,T,L^\infty(\mathbb{R}^N))} \right) < \infty, $$

$$ \bar{A}_1 := \|\nabla_1 K\|_{L^\infty(\mathbb{R}^N \times \mathbb{R}^N)} \|\partial_1 C(d,p)\|_{L^\infty(0,T,L^\infty(\mathbb{R}^N))} \|\eta\|_{L^\infty(0,T,L^\infty(\mathbb{R}^N))} \exp \left( \|g(d,p)\|_{L^1(0,T,L^\infty(\mathbb{R}^N))} \right) < \infty, $$

$$ A_2 := \|p\|_{C(0,T;L^1(\mathbb{R}^N))} \|\operatorname{div}(V')\|_{C(0,T;L^\infty(\mathbb{R}^N))} + 2 \|V'\|_{C(0,T;L^\infty(\mathbb{R}^N))} \|\nabla p\|_{C(0,T;L^1(\mathbb{R}^N))} < \infty, $$

are positive constants. We get after some usual estimates:

$$ \left| L(w)(t, x) \right| \le T \left( A_1 + A_2 \|K\|_{L^\infty(\mathbb{R}^N \times \mathbb{R}^N)} \exp \left( \|g(d,p)\|_{L^1(0,T,L^\infty(\mathbb{R}^N))} \right) \|w\|_{C(0,T;C_b^1(\mathbb{R}^N))} \right), $$

and

$$ \left| \nabla L(w)(t, x) \right| \le T \left( \bar{A}_1 + A_2 \|\nabla K\|_{L^\infty(\mathbb{R}^N \times \mathbb{R}^N)} \exp \left( \|g(d,p)\|_{L^1(0,T,L^\infty(\mathbb{R}^N))} \right) \|w\|_{C(0,T;C_b^1(\mathbb{R}^N))} \right). $$
Therefore, by taking the supremum over time and space, and choosing $T$ small enough we get:

$$ \|\mathrm{L}(w)\|_{C(0,T;C_b^1(\mathbb{R}^N))} \leq M. \qquad (40) $$

Hence, we deduce that the mapping $\mathrm{L}$ is well-defined. Now let $w_1, w_2 \in B$, we get using similar estimations:

$$ \|\mathrm{L}(w_1) - \mathrm{L}(w_2)\|_{C(0,T;C_b^1(\mathbb{R}^N))} \leq TR \|w_1 - w_2\|_{C(0,T;C_b^1(\mathbb{R}^N))}. \qquad (41) $$

where $R$ is a constant depending on $K, g, V$ and $\eta$. Hence, for $T$ small enough, we deduce that $L$ is a contacton on $B$, which is a closed subset of the Banach space $C(0, T; C_b^1(\mathbb{R}^N))$, then by Banach's fixed point theorem, there exists a unique function $w^* \in B$ such that $w^* = \mathrm{L}[w^*]$. Substituting back into the equation (36), and by following the approach of [11], we get the existence of a solution for every $T > 0$.

For the uniqueness, let $(\xi_1, \eta_1)$ and $(\xi_2, \eta_2)$ be two solutions of the linearized system. Let define $\tilde{\xi} = \xi_1 - \xi_2$, and $\tilde{\eta} = \eta_1 - \eta_2$. Then for the parabolic equation $\tilde{\eta}$ we have:

$$ \begin{cases} \frac{\partial \tilde{\eta}}{\partial t} - \mathrm{div}(D(x)\nabla \tilde{\eta}) + \partial_2 \Gamma(I, d) \tilde{\eta} = 0, \\ \tilde{\eta}(0, x) = 0. \end{cases} $$

By using estimates from (3.1), we deduce that $\tilde{\eta} = 0$ in $C(0, T, H^1(\mathbb{R}^N))$. In the other hand, for $\tilde{\xi}$ we have:

$$ \begin{cases} \frac{\partial \tilde{\xi}}{\partial t} + \mathrm{div} \left( V[w_p] \tilde{\xi} \right) + \mathrm{div} \left( V'[w_p] w_{\tilde{\xi}} p \right) = F'(p) \tilde{\xi} - \partial_2 C(d, p) \tilde{\xi}, \\ \tilde{\xi}(0, x) = 0. \end{cases} $$

By using the characteristics as in (37), we obtain:

$$ \tilde{\xi}(t, x) = - \int_{0}^{t} \mathrm{div} \left( V'[w_p] w_{\tilde{\xi}} p \right) (s, X[t, x](s)) \exp \left( \int_{s}^{t} g(d, p)(\tau, X[t, x](\tau)) d\tau \right) \det(D_x X[t, x](s)) ds, $$

By taking the $L^1$ estimate as before we get:

$$ \|\tilde{\xi}(t, \cdot)\|_{L^1(\mathbb{R}^N)} \leq \int_0^t \|\mathrm{div}(V'[w_p]w_{\tilde{\xi}}p)\|_{L^1(\mathbb{R}^N)} \exp \left( \|g(d,p)\|_{L^\infty(\mathbb{R}^N)}(t-s) \right) ds. $$

Using assumptions (A8) and (A9):

$$ \|\mathrm{div}(V'[w_p]w_{\tilde{\xi}}p)\|_{L^1(\mathbb{R}^N)} \leq C_1 \|\tilde{\xi}(t, \cdot)\|_{L^1(\mathbb{R}^N)}, \qquad (42) $$

where $C_1$ depends on $\|K\|_{C_b^1(\mathbb{R}^N)}$, $\|V'\|_{C_b^1(\mathbb{R}^N)}$, and $\|p\|_{L^\infty(0,T,L^\infty(\mathbb{R}^N))}$. Therefore:

$$ \|\tilde{\xi}(t, \cdot)\|_{L^1(\mathbb{R}^N)} \leq C_1 \int_0^t \|\tilde{\xi}(s, \cdot)\|_{L^1(\mathbb{R}^N)} \exp(C_2(t-s)) ds. $$

By Gronwall's inequality we have $\|\tilde{\xi}(t, \cdot)\|_{L^1(\mathbb{R}^N)} \leq 0$, which leads to the desired result. $\square$
**Lemma 4.2** (Gâteaux Differentiability of the Control-to-State Mapping). Let Assumptions (A1)-(A10) hold. Then G given by (34) is Gâteaux differentiable at any $I \in \mathcal{U}_{ad}$. Furthermore, for any direction $h \in L^\infty(0,T)$, the Gâteaux derivative of G at I in the direction h is $\frac{\partial G}{\partial h}(I) = (\xi, \eta)$, where $(\xi, \eta)$ is the solution to the linearized system (35).

*Proof.* Let $I \in \mathcal{U}_{ad}$, $h \in L^\infty(0,T)$, and $\varepsilon > 0$ such that $I^\varepsilon := I + \varepsilon h \in U_{ad}$. Let $(p^\varepsilon, d^\varepsilon) = G(I^\varepsilon)$ and $(p, d) = G(I)$, and let define $\xi^\varepsilon = (p^\varepsilon - p)/\varepsilon$ and $\eta^\varepsilon = (d^\varepsilon - d)/\varepsilon$. We aim to show that $(\xi^\varepsilon, \eta^\varepsilon) \to (\xi, \eta)$ in X as $\varepsilon \to 0$. We have for difference quotients $(\xi^\varepsilon, \eta^\varepsilon)$:

$$ 
\begin{aligned}
\frac{\partial \xi^\varepsilon}{\partial t} + \operatorname{div} \left( \frac{V[w_{p^\varepsilon}] p^\varepsilon - V[w_p] p}{\varepsilon} \right) &= \frac{F(p^\varepsilon) - F(p)}{\varepsilon} - \frac{C(d^\varepsilon, p^\varepsilon) - C(d, p)}{\varepsilon}, \\
\frac{\partial \eta^\varepsilon}{\partial t} - \operatorname{div}(D(x)\nabla \eta^\varepsilon) &= \frac{\Gamma(I^\varepsilon, d^\varepsilon) - \Gamma(I, d)}{\varepsilon}, \\
\xi^\varepsilon(x, 0) = 0, \quad \eta^\varepsilon(x, 0) = 0.
\end{aligned}
$$

We have $w_{p^\varepsilon} = w_p + \varepsilon w_{\xi^\varepsilon}$. By assumptions (A9) and (A10), Taylor's theorem with integral remainder for the nonlinear terms gives:

$$ 
\begin{aligned}
\frac{V[w_{p^\varepsilon}] p^\varepsilon - V[w_p] p}{\varepsilon} &= V[w_p] \xi^\varepsilon + V'[w_p] w_{\xi^\varepsilon} p + R_1^\varepsilon, \\
\frac{F(p^\varepsilon) - F(p)}{\varepsilon} &= F'(p)\xi^\varepsilon + R_2^\varepsilon, \\
\frac{C(d^\varepsilon, p^\varepsilon) - C(d, p)}{\varepsilon} &= \partial_1 C(d, p)\eta^\varepsilon + \partial_2 C(d, p)\xi^\varepsilon + R_3^\varepsilon, \\
\frac{\Gamma(I^\varepsilon, d^\varepsilon) - \Gamma(I, d)}{\varepsilon} &= \partial_1 \Gamma(I, d)h + \partial_2 \Gamma(I, d)\eta^\varepsilon + R_4^\varepsilon,
\end{aligned}
$$

where $\|R_i^\varepsilon\|_{L^\infty(0,T;L^\infty(\mathbb{R}^N))} = o(\varepsilon)$, as $\varepsilon \to 0$, for $i = 1, \dots, 4$. Therefore we have:

$$ 
\left\{
\begin{aligned}
\frac{\partial \xi^\varepsilon}{\partial t} + \operatorname{div}(V[w_p] \xi^\varepsilon) + \operatorname{div}(V'[w_p] w_{\xi^\varepsilon} p) &= F'(p)\xi^\varepsilon - \partial_1 C(d, p)\eta^\varepsilon - \partial_2 C(d, p)\xi^\varepsilon - \operatorname{div}(R_1^\varepsilon) + R_2^\varepsilon - R_3^\varepsilon, && (a) \\
\frac{\partial \eta^\varepsilon}{\partial t} - \operatorname{div}(D(x)\nabla \eta^\varepsilon) &= \partial_1 \Gamma(I, d)h + \partial_2 \Gamma(I, d)\eta^\varepsilon + R_4^\varepsilon, && (b) \\
\xi^\varepsilon(0, x) = 0, \quad \eta^\varepsilon(0, x) = 0.
\end{aligned}
\right.
\tag{43}
$$

Now let $\tilde{\xi} = \xi^\varepsilon - \xi$ and $\tilde{\eta} = \eta^\varepsilon - \eta$. By subtracting (35) from (43) we get the following perturbed system:

$$ 
\left\{
\begin{aligned}
\frac{\partial \tilde{\xi}}{\partial t} + \operatorname{div}(V[w_p] \tilde{\xi}) + \operatorname{div}(V'[w_p] p w_{\tilde{\xi}}) &= g(d, p) \tilde{\xi} - \partial_1 C(d, p) \tilde{\eta} - \operatorname{div}(R_1^\varepsilon) + R_2^\varepsilon - R_3^\varepsilon, && (a) \\
\frac{\partial \tilde{\eta}}{\partial t} - \operatorname{div}(D(x) \nabla \tilde{\eta}) &= \partial_2 \Gamma(I, d) \tilde{\eta} + R_4^\varepsilon, && (b) \\
\tilde{\xi}(x, 0) = 0, \quad \tilde{\eta}(x, 0) = 0.
\end{aligned}
\right.
$$

By using a similar technique as in the uniqueness proof of linearized system, we obtain:

$$ \|(\tilde{\xi}, \tilde{\eta})\|_X \le C_1 \|R^\varepsilon\|_{L^\infty(0,T;L^\infty(\mathbb{R}^N))}, $$
where $R^\epsilon = R_2^\epsilon + \text{div}(R_1^\epsilon) + R_3^\epsilon + R_4^\epsilon$ and $C_1$ is some positive constant independent of $\epsilon$. Therefore,

$$ \lim_{\epsilon \to 0} \left\| \frac{G(I + \epsilon h) - G(I)}{\epsilon} - (\xi, \eta) \right\|_X = \lim_{\epsilon \to 0} \left\| (\tilde{\xi}, \tilde{\eta}) \right\|_X = 0. $$

This establishes the Gâteaux differentiability of $G$ at $I$ in the direction $h$, with $G'(I)h = (\xi, \eta)$. $\square$

Now, we establish a result on the differentiability of the cost functional:

**Theorem 4.1** (On the Gâteau-differentiability of the Cost Functional). Let Assumptions (A1)-(A10) hold. Then the cost functional $J$ is Gâteaux differentiable at any $I \in \mathcal{U}_{ad}$. Furthermore, for any direction $h \in L^\infty(0,T)$, the Gâteaux derivative of $J$ at $I$ in the direction $h$ is given by:

$$ \langle J'(I), h \rangle = \int_0^T \left[ \alpha \int_{\mathbb{R}^N} \xi(t,x)dx + 2\beta I(t)h(t) \right] dt + \gamma \int_{\mathbb{R}^N} \xi(T,x)dx \quad (45) $$

where $(\xi, \eta)$ is the solution to the linearized system (35).

*Proof.* Let $I \in \mathcal{U}_{ad}$ and define the functional $L_I : L^\infty(0,T) \to \mathbb{R}$ by:

$$ \langle L_I, h \rangle = \int_0^T \left[ \alpha \int_{\mathbb{R}^N} \xi(t,x)dx + 2\beta I(t)h(t) \right] dt + \gamma \int_{\mathbb{R}^N} \xi(T,x)dx \quad (46) $$

where $(\xi, \eta)$ solves the linearized system (35) with the direction $h$. Now let $h_1, h_2 \in L^\infty(0,T)$ and $\lambda \in \mathbb{R}$, and let us consider $h = h_1 + \lambda h_2$, with $(\xi_1, \eta_1)$ and $(\xi_2, \eta_2)$ are the solutions to the linearized system corresponding to directions $h_1$ and $h_2$ respectively. By the linearity of system (35), the solution $(\xi, \eta)$ corresponding to $h$ satisfies $\xi = \xi_1 + \lambda \xi_2$ and $\eta = \eta_1 + \lambda \eta_2$. Consequently:

$$ \langle L_I, h_1 + \lambda h_2 \rangle = \int_0^T \left[ \alpha \int_{\mathbb{R}^N} (\xi_1 + \lambda \xi_2)dx + 2\beta I(t)(h_1(t) + \lambda h_2(t)) \right] dt + \gamma \int_{\mathbb{R}^N} (\xi_1 + \lambda \xi_2)(T,x)dx \\ = \langle L_I, h_1 \rangle + \lambda \langle L_I, h_2 \rangle. $$

Therefore $L_I$ is linear. Furthermore, for any $h \in L^\infty(0,T)$, by using the estimates from Lemma 4.1 we obtain:

$$ |\langle L_I, h \rangle| \le \alpha \int_0^T \int_{\mathbb{R}^N} |\xi(t,x)|dxdt + 2\beta \|I\|_{L^\infty} \|h\|_{L^\infty} T + \gamma \int_{\mathbb{R}^N} |\xi(T,x)|dx \le C_1 \|h\|_{L^\infty(0,T)}, $$

where $C_1$ is a positive constant depending on $\alpha, \beta, \gamma, T$, and the bounds from Lemma 4.1. Therefore $L_I$ is continuous.

To show that $L_I$ represents the Gâteaux derivative of $J$, let $h \in L^\infty(0,T)$, $\epsilon > 0$, and let define $I^\epsilon = I + \epsilon h$ and let $(p^\epsilon, d^\epsilon) = G(I^\epsilon)$ and $(p, d) = G(I)$. Let define $\Phi_I(\epsilon) = \frac{J(I^\epsilon) - J(I)}{\epsilon}$, we have:

$$ \Phi_I(\epsilon) = \int_0^T \left[ \alpha \int_{\mathbb{R}^N} \frac{p^\epsilon - p}{\epsilon} dx + \beta \frac{(I^\epsilon)^2 - I^2}{\epsilon} \right] dt + \gamma \int_{\mathbb{R}^N} \frac{p^\epsilon(T,x) - p(T,x)}{\epsilon} dx \\ = \int_0^T \left[ \alpha \int_{\mathbb{R}^N} \xi^\epsilon dx + \beta(2Ih + \epsilon h^2) \right] dt + \gamma \int_{\mathbb{R}^N} \xi^\epsilon(T,x) dx $$

where $\xi^\epsilon = (p^\epsilon - p)/\epsilon$. By Lemma 4.1, we have $\xi^\epsilon \to \xi$ in $X$ as $\epsilon \to 0$, where $\xi$ is the solution to the linearized system (35). Therefore:

$$ \lim_{\epsilon \to 0} \Phi_I(\epsilon) = \int_0^T \left[ \alpha \int_{\mathbb{R}^N} \xi(t,x)dx + 2\beta I(t)h(t) \right] dt + \gamma \int_{\mathbb{R}^N} \xi(T,x)dx = \langle L_I, h \rangle. \quad (47) $$

This establishes that $J$ is Gâteaux differentiable at $I$ with derivative $J'(I) = L_I$, which complete the proof. $\square$
## 4.2. Necessary conditions

To simplify the expression of the gradient of $J$, we introduce an adjoint system. This system is designed to eliminate the dependence on $\xi$ in the expression of $J'(I)$. We introduce the following system:

$$ egin{cases} 
-\frac{\partial q}{\partial t}(t, x) - V[w_p(t, x)] \cdot \nabla q(t, x) - \int_{\Omega} K(y, x) (pV'[w_p] \cdot \nabla q)(t, y) dy + g(d, p)(t, x)q(t, x) = \alpha, & (a) \\ 
-\frac{\partial r}{\partial t}(t, x) - \operatorname{div}(D(x)\nabla r(t, x)) + \partial_2 \Gamma(I, d)(t, x)r(t, x) = -\partial_2 C(d, p)q(t, x), & (b) \\ 
q(T, x) = \gamma, \quad r(T, x) = 0. 
\end{cases} 
\quad (48) $$

We first establish the following lemma:

**Lemma 4.3** (Well-posedness of the Adjoint System). Let Assumptions (A1)-(A12) hold. Then the adjoint system (48) has a unique solution $(q, r)$ in $C(0, T; L^{\infty}(\mathbb{R}^N)) \times C(0, T; H^1(\mathbb{R}^N))$.

*Proof.* Let $s = T - t$ and define:

$$ \bar{q}(s, x) = q(T - s, x), \quad \bar{r}(s, x) = r(T - s, x). \quad (49) $$

The transformed system becomes:

$$ egin{cases} 
\frac{\partial \bar{q}}{\partial s} - V[w_p] \cdot \nabla \bar{q} - \int_{\Omega} K(y, x) (pV'[w_p] \cdot \nabla \bar{q})(T-s, y) dy + g(d, p)\bar{q} = \alpha, & (a) \\ 
\frac{\partial \bar{r}}{\partial s} - \operatorname{div}(D(x)\nabla \bar{r}) + \partial_2 \Gamma(I, d)\bar{r} = -\partial_2 C(d, p)\bar{q}, & (b) \\ 
\bar{q}(0, x) = \gamma, \quad \bar{r}(0, x) = 0. 
\end{cases} 
\quad (50) $$

By using assumptions (A11) and (A12), the integration by parts in the nonlocal term leads to:

$$ \int_{\Omega} K(y, x) (pV'[w_p] \cdot \nabla \bar{q})(T-s, y) dy = - \int_{\Omega} \bar{q}(s, y) \operatorname{div}_y (p(y)V'[w_p(y)]K(y, x)) dy, \quad (51) $$

where the boundary terms vanish since $K(y, x) = 0$ for $y \in \partial\Omega$ by (A10) and $V'(0) = 0$ by (A11). Therefore, $(\bar{q}, \bar{r})$ satisfies:

$$ egin{cases} 
\frac{\partial \bar{q}}{\partial s} - \operatorname{div}(V[w_p]\bar{q}) + \int_{\Omega} \bar{q}(y) \operatorname{div}_y (p(y)V'[w_p(y)]K(y, x)) dy = -(g(d, p) + \operatorname{div}(V[w_p])) \bar{q} + \alpha, \\ 
\frac{\partial \bar{r}}{\partial s} - \operatorname{div}(D(x)\nabla \bar{r}) + \partial_2 \Gamma(I, d)\bar{r} = -\partial_2 C(d, p)\bar{q}, \\ 
\bar{q}(0, x) = \gamma, \quad \bar{r}(0, x) = 0. 
\end{cases} 
\quad (52) $$

The transformed system has a similar structure as the linearized system (35). Therefore, by the same argument as in Lemma 4.1, the system (52) has a unique solution $(\bar{q}, \bar{r})$, which imply the existence and uniqueness of $(q, r)$ for the adjoint system (48). □

Now we can state and prove the necessary optimality conditions:
**Theorem 4.2** (Necessary Optimality Conditions). Let $I^* \in \mathcal{U}_{ad}$ be an optimal control for the problem with corresponding state $(p^*, d^*)$. Then there exists an adjoint state $(q, r)$ satisfying the adjoint system (48) such that:

$$ \int_{0}^{T} \left( 2\beta I^{*}(t) + \int_{\mathbb{R}^{N}} \partial_{1}\Gamma(I^{*}, d^{*})r(t, x)dx \right) (v(t) - I^{*}(t)) dt \geq 0 \quad (53) $$

for all $v \in \mathcal{U}_{ad}$.

*Proof.* Let $I^*$ be an optimal control. For any $v \in \mathcal{U}_{ad}$ and $\varepsilon \in (0, 1)$, define $I^\varepsilon = I^* + \varepsilon(v - I^*)$. Note that $I^\varepsilon \in \mathcal{U}_{ad}$ due to the convexity of $\mathcal{U}_{ad}$. By the optimality of $I^*$, we have:

$$ 0 \le \lim_{\varepsilon \to 0^+} \frac{J(I^\varepsilon) - J(I^*)}{\varepsilon} = \langle J'(I^*), v - I^* \rangle. \quad (54) $$

From Theorem 4.1, we have:

$$ \langle J'(I^*), v - I^* \rangle = \int_{0}^{T} \left[ \alpha \int_{\mathbb{R}^{N}} \xi(t, x)dx + 2\beta I^{*}(t)(v(t) - I^{*}(t)) \right] dt + \gamma \int_{\mathbb{R}^{N}} \xi(T, x)dx \quad (55) $$

where $(\xi, \eta)$ solves the linearized system with $h = v - I^*$.

We multiply the first equation of the linearized system (35-a) by $q$, and the second equation (35-b) by $r$, then integrate over space and time:

$$ \int_{0}^{T} \int_{\mathbb{R}^{N}} q \left( \frac{\partial \xi}{\partial t} + \operatorname{div}(V[w_{p^*}]\xi) + \operatorname{div}(V'[w_{p^*}](w_{\xi})p^*) - g(d^*, p^*)\xi + \partial_1 C(d^*, p^*)\eta \right) dxdt = 0, $$

$$ \int_{0}^{T} \int_{\mathbb{R}^{N}} r \left( \frac{\partial \eta}{\partial t} - \operatorname{div}(D(x)\nabla \eta) - \partial_{2}\Gamma(I, d)\eta - \partial_{1}\Gamma(I, d)(v - I^{*}) \right) dxdt = 0. $$

After integration by parts we get for $\xi$ equation:

$$ \begin{aligned} & \gamma \int_{\mathbb{R}^N} \xi(T, x)dx - \int_0^T \int_{\mathbb{R}^N} \frac{\partial q}{\partial t} \xi dxdt - \int_0^T \int_{\mathbb{R}^N} \nabla q \cdot (V[w_{p^*}]\xi) dxdt \\ & - \int_0^T \int_{\mathbb{R}^N} \nabla q \cdot (V'[w_{p^*}]w_{\xi}p^*) dxdt \\ & = \int_0^T \int_{\mathbb{R}^N} qF'(p^*)\xi dxdt - \int_0^T \int_{\mathbb{R}^N} q\partial_2 C\xi dxdt - \int_0^T \int_{\mathbb{R}^N} q\partial_1 C\eta dxdt. \end{aligned} $$

We have by using Fubini theorem:

$$ \begin{aligned} \int_{0}^{T} \int_{\mathbb{R}^{N}} \nabla q(t, x) \cdot (V'[w_{p^*}(t, x)]w_{\xi}(t, x)p^*(t, x)) dxdt &= \int_{0}^{T} \int_{\mathbb{R}^{N}} \int_{\Omega} \nabla q(t, x) \cdot V'[w_{p^*}(t, x)]k(x, y)\xi(t, y)dxdydt \\ &= \int_{0}^{T} \int_{\mathbb{R}^{N}} \xi(t, x) \int_{\Omega} \nabla q(t, y) \cdot (V'[w_{p^*}(t, y)]k(y, x))dydxdt. \end{aligned} $$

Hence we get:

$$ \begin{aligned} & \gamma \int_{\mathbb{R}^N} \xi(T, x)dx + \int_0^T \int_{\mathbb{R}^N} \xi \left[ -\frac{\partial q}{\partial t} - \nabla q \cdot (V[w_{p^*}]) - \int_{\mathbb{R}^N} \nabla q \cdot (V'[w_{p^*}]p^*k(x, \cdot))dx - g(d^*, p^*)q \right] dxdt \\ & - \int_0^T \int_{\mathbb{R}^N} q\partial_1 C(d^*, p^*)\eta dxdt = 0. \end{aligned} $$
By substituting (48) we get:

$$ \gamma \int_{\mathbb{R}^N} \xi(T, x)dx + \alpha \int_0^T \int_{\mathbb{R}^N} \xi(t, x)dxdt + \int_0^T \int_{\mathbb{R}^N} q\partial_1 C(d^*, p^*)\eta dxdt = 0. \quad (56) $$

In the other hand, we have for the second equation after integrating by parts:

$$ \begin{aligned} & \int_{\mathbb{R}^N} r(T, x)\eta(T, x)dx - \int_{\mathbb{R}^N} r(x, 0)\eta(x, 0)dx - \int_0^T \int_{\mathbb{R}^N} \frac{\partial r}{\partial t}\eta dxdt \\ & + \int_0^T \int_{\mathbb{R}^N} \nabla r \cdot (D(x)\nabla\eta)dxdt - \int_0^T \int_{\mathbb{R}^N} r\partial_2\Gamma(I^*, d^*)\eta dxdt \\ & = \int_0^T \int_{\mathbb{R}^N} r\partial_1\Gamma(I^*, d^*)(v - I^*)dxdt, \end{aligned} $$

after substituting (48) we get:

$$ \int_0^T \int_{\mathbb{R}^N} q\partial_1 C(d^*, p^*)\eta dxdt = - \int_0^T \int_{\mathbb{R}^N} r\partial_1 \Gamma(I^*, d^*)(v - I^*)dxdt. \quad (57) $$

By substituting (56) and (57) back into (45) we obtain:

$$ \int_0^T \left( 2\beta I^*(t) + \int_{\mathbb{R}^N} \partial_1 \Gamma(I^*, d^*) r(x, t) dx \right) (v(t) - I^*(t)) dt \ge 0 \quad (58) $$

for all $v \in \mathcal{U}_{ad}$, which completes the proof. $\square$

**Corollary 4.1** (Characterization of Optimal Control). Let Assumptions (A1)-(A12) hold and let $I^* \in \mathcal{U}_{ad}$ be an optimal control with corresponding adjoint state $r$. Then:

$$ I^*(t) = \max \left\{ 0, \min \left\{ M, -\frac{1}{2\beta} \int_{\mathbb{R}^N} \partial_1 \Gamma(I^*, d^*) r(t, x) dx \right\} \right\} \quad \text{a.e. } t \in [0, T]. \quad (59) $$

*Proof.* Let $I^*$ be an optimal control for problem ((5)). From Theorem 4.2, we have the following variational inequality:

$$ \int_0^T \left( 2\beta I^*(t) + \int_{\mathbb{R}^N} \partial_1 \Gamma(I^*, d^*) r(t, x) dx \right) (v(t) - I^*(t)) dt \ge 0, \quad \forall v \in \mathcal{U}_{ad}. \quad (60) $$

We introduce the switching function $\psi(t) := 2\beta I^*(t) + \int_{\mathbb{R}^N} \partial_1 \Gamma(I^*, d^*) r(t, x) dx$. The variational inequality can thus be rewritten as:

$$ \int_0^T \psi(t)(v(t) - I^*(t))dt \ge 0, \quad \forall v \in \mathcal{U}_{ad}. \quad (61) $$

Let us consider first the case where $\psi(t) > 0$ for some $t \in [0, T]$. Let suppose, by contradiction, that there exists a set of positive measure where $\psi(t) > 0$ and $I^*(t) > 0$. For sufficiently small $\varepsilon > 0$, we can define an admissible control $v_\varepsilon(t) := \max\{0, I^*(t) - \varepsilon\}$.

For this choice of $v_\varepsilon$, we have $v_\varepsilon(t) - I^*(t) = -\varepsilon < 0$, which leads to $\int_0^T \psi(t)(v_\varepsilon(t) - I^*(t))dt < 0$, which contradicts the optimality condition, thus $I^*(t) = 0$ whenever $\psi(t) > 0$.
For the case where $\psi(t) < 0$, we can apply a similar argument. Let suppose there exists a set of positive measure where $\psi(t) < 0$ and $I^*(t) < M_{tol}$. We define $v_{\epsilon}(t) = \min\{M_{tol}, I^*(t) + \epsilon\}$. This choice leads to a similar contradiction, establishing that $I^*(t) = M_{tol}$ whenever $\psi(t) < 0$.

When $\psi(t) = 0$, we have directly:

$$2\beta I^*(t) + \int_{\mathbb{R}^N} \partial_1 \Gamma(I^*, d^*) r(t, x) dx = 0.$$

Therefore:

$$I^*(t) = -\frac{1}{2\beta} \int_{\mathbb{R}^N} \partial_1 \Gamma(I^*, d^*) r(t, x) dx.$$

The constraints $0 \le I^*(t) \le M_{tol}$ must hold for all $t \in [0, T]$. Combining these constraints with our previous analysis, we obtain the complete characterization:

$$I^*(t) = \max \left\{ 0, \min \left\{ M_{tol}, -\frac{1}{2\beta} \int_{\mathbb{R}^N} \partial_1 \Gamma(I^*, d^*) r(t, x) dx \right\} \right\} \quad \text{a.e. } t \in [0, T]. \qquad (62)$$

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

# 5. Numerical Simulations

For numerical implementation of the optimal control strategy characterized by equation (59), we formulate a computational approach based on the variational structure of the control problem and the adjoint system derived in Section 4. The iterative scheme is as follows:

* At iteration $k$, given the current control approximation $I^k \in U_{ad}$.

- Compute the state variables $(p^k, d^k)$ by solving the coupled system (6).

- Compute the adjoint variables $(q^k, r^k)$ through backward integration of system (48).

- Evaluate the gradient of the cost functional $\nabla J(I^k)$.

* Update $I^{k+1} = \Pi_{U_{ad}} (I^k - \alpha_k \nabla J(I^k))$, where $\Pi_{U_{ad}}$ denotes the projection operator onto the admissible set $U_{ad}$ and $\alpha_k > 0$ represents the step size.

While this iterative procedure provides a practical realization of the theoretical framework established in preceding sections, a formal convergence analysis of the numerical scheme lies beyond the scope of the present work. The numerical simulations were implemented using Python 3.13 with NumPy library.

## 5.1. Discretization and Parameter Configuration

The computational domain $\Omega = [0, 1]$ was discretized with a uniform spatial grid of $N_x = 200$ points, while the time domain $[0, T]$ with $T = 1$ was discretized according to the Courant-Friedrichs-Lewy (CFL) condition with $CFL = 0.025$. Let $\{x_i\}_{i=1}^{N_x}$ and $\{t_n\}_{n=1}^{N_t}$ denote the spatial and temporal discretization points, respectively. The initial condition $p_0(x)$ was specified as a normalized Gaussian function centered at $x = 0.5$ with standard deviation 0.1, while $d_0(x) \equiv 0$ throughout the domain. For the abstract functions introduced in system (6), we specify the nonlocal velocity field $V[w_p] = \kappa \cdot w_p(t, x)$, the logistic growth function $F(p) = r \cdot p \cdot (1 - p/K_p)$, the drug-induced cell death term $C(d, p) = \delta \cdot d \cdot p$, and the drug exchange and clearance mechanism
$\Gamma(I, d) = \Gamma \cdot (I - d) - \lambda \cdot d$ (see [19]). Where $w_p(t, x) = \int_{\Omega} K(x, y)p(t, y) dy$ represents the nonlocal term with $K$ defined as a Gaussian kernel with standard deviation $\sigma$.

For the spatial discretization, the advection term in the equation (6)-(b) was approximated using a Lax-Friedrichs scheme, the diffusion term in the equation (6)-(c) was discretized using a standard second-order central difference approximation, and the nonlocal convolution integral (3) was computed using the Fast Fourier Transform (FFT). The resulting semi-discrete system was integrated in time using an explicit Euler method with time step constrained by the CFL condition.

## 5.2. Results and Analysis

The parameter values used in simulations were: $r = 1.2$, $K_p = 1.0$, $\delta = 0.8$, $\kappa = 0.05$, $\sigma = 0.01$, $D = 0.05$, $\Gamma = 2.5$, and $\lambda = 0.3$. For the cost functional defined in equation (4), the weights were $\alpha = 1.0$, $\beta = 0.1$, and $\gamma = 1.0$, with the maximum tolerable drug concentration constrained to $M_{\text{tol}} = 4.0$.

<figure><img src="image_3.png" /><figcaption>**Figure 1:** Optimal control profile $I^*(t)$ over the treatment interval $[0, T]$.</figcaption></figure>

The computed optimal control $I^*(t)$ is a monotonically decreasing function over the treatment interval, as shown in Figure 1. This decay is pronounced during the initial phase $t \in [0, 0.3]$, followed by a more gradual reduction, suggesting that the optimal strategy entails substantial initial drug administration followed by progressive tapering. The asymptotic approach of $I^*(t)$ toward zero as $t \to T$ reflects the balance between tumor reduction efficacy and control cost in the objective functional.

The evolution of tumor density under optimal control shows a differences compared to the uncontrolled case, as illustrated in Figure 2 and 4. Let $p_c(t, x)$ and $p_u(t, x)$ denote the tumor density functions under controlled and uncontrolled conditions, respectively. The relation $\sup_{x \in \Omega} p_c(t, x) < \sup_{x \in \Omega} p_u(t, x)$ holds for all $t \in (0, T]$. The spatial profiles at discrete time points (Figure 3) reveal progressive divergence between the controlled and uncontrolled distributions, having in approximately 27% reduction in peak tumor density by $t = 0.9$.
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2:** Three-dimensional visualization of tumor density evolution $p(t, x)$.</figcaption></figure>

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

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

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

<figure><img src="image_9.png" /><figcaption>**Figure 3:** Tumor density profiles $p(t, x)$ at selected time points comparing controlled (solid blue) and uncontrolled (dashed red) cases.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 4: Contour plots of tumor density $p(t, x)$ in the $(x, t)$ plane.</figcaption></figure>

The distribution of drug concentration $d(t, x)$ (Figure 5) shows that this concentration is reaching a maximum value at $t \approx T/2$ despite the monotonically decreasing control input. This phase lag phenomenon emerges from the interaction between the source term and the diffusion operator, with the spatial homogeneity of $d(t, x)$ in contrast to the heterogeneous distribution of tumor density.

<figure><img src="image_11.png" /><figcaption>Figure 5: Three-dimensional visualization of drug concentration $d(t, x)$.</figcaption></figure>

To quantify the efficacy of the optimal control strategy, we define functionals $M_c(t) = \int_{\Omega} p_c(t, x) dx$ and $M_u(t) = \int_{\Omega} p_u(t, x) dx$, representing the tumor mass under optimal control and without control, respectively. The relative improvement metric $\Delta(t) = \frac{M_u(t)-M_c(t)}{M_u(t)} \times 100\%$ shows an increasing behavior for $t \in [0, T]$, with $\lim_{t \to T} \Delta(t) \approx 27\%$ (Figure 6), which quantifies the reduction achieved through the implementation of optimal control.
<figure><img src="image_12.png" /><figcaption>**Figure 6:** Relative improvement $\Delta(t)$ in tumor reduction due to optimal control intervention over time.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7:** Comparison of tumor mass evolution with optimal control (solid blue) and without control (dashed red).</figcaption></figure>

The numerical optimization procedure shows stable convergence characteristics, shown by the monotonic decrease in both the cost functional $J(I^k)$ and the gradient norm $\|\nabla J(I^k)\|_{L^2(0,T)}$ (Figure 8). The algorithm was executed for a maximum of 50 iterations with the convergence criteria $|J(I^{k+1}) - J(I^k)|/|J(I^k)| < 5 \times 10^{-5}$.

<figure><img src="image_14.png" /><figcaption>**Figure 8:** Evolution of optimization metrics across iterations: (a) monotonic decrease in cost functional $J(I^k)$ and (b) corresponding reduction in gradient norm $\|\nabla J(I^k)\|_{L^2(0,T)}$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 9:** Evolution of control profile from initial guess (dashed gray) to optimal solution (solid black).</figcaption></figure>

The initial control approximation, an exponential decay profile $I^0(t) = M_{tol}e^{-4t}$, was chosen to align with the decreasing structure predicted by the necessary optimality conditions in equation (59) while starting at the maximum tolerable dose. As illustrated in Figure 9, this initial guess was progressively refined through optimization, converging to a profile with a steeper initial decrease and a smoother decrease.

## Conclusion

The development of effective cancer treatment strategies remains a critical challenge in modern medicine, requiring sophisticated mathematical frameworks that can capture the complex dynamics of tumor growth while providing practical guidance for therapy optimization. In this work, we have presented a mathematical analysis of a nonlocal tumor growth model incorporating chemotherapy, making some contributions to the theoretical understanding of optimal drug delivery.

Our analysis began by establishing a theoretical framework for studying nonlocal tumor growth models with chemotherapy, proving the existence and uniqueness of a solution under some assumptions. The incorporation of nonlocal terms, representing long-range cellular interactions, provides a more realistic model of tumor dynamics than traditional local approaches. We then developed an optimal control theory for this system, deriving the optimality conditions through a novel adjoint formulation that accounts for the nonlocal effects.

We obtained an explicit characterization of the optimal control strategy that provides practical guidance for drug dosing protocols. Our analysis reveals that optimal treatment schedules follow a precise structure determined by the interplay between tumor dynamics and drug effects. This structure, characterized by a switching function involving the adjoint state, suggests that effective chemotherapy regimens should adapt to both the current tumor state and its predicted evolution.

The mathematical techniques developed in this work extend beyond oncology, with potential applications in numerous biological systems where nonlocal interactions govern emergent behaviors. These include neural tissue dynamics, bacterial biofilm formation, embryogenesis, and even ecosystem population distributions where long-range interactions shape spatial patterns. Within cancer research specifically, several promising extensions warrant investigation. Future work could incorporate immune system responses by coupling our model with immune cell dynamics, allowing exploration of combined chemo-immunotherapy protocols. The framework could also be extended to account for tumor heterogeneity and phenotypic plasticity, including drug resistance mechanisms that evolve during treatment.
Acknowledgment

This work was supported by the Centre National pour la Recherche Scientifique et Technique (CNRST) of Morocco and the Laboratoire de Mathématiques Pures et Appliquées (LMPA) at the Université du Littoral Côte d'Opale (ULCO), France.

Bibliography

[1] P. M. Altrock, L. L. Liu, and F. Michor. The mathematics of cancer: integrating quantitative models. *Nature Reviews Cancer*, 15(12):730-745, 2015.

[2] N. J. Armstrong, K. J. Painter, and J. A. Sherratt. A continuum approach to modelling cell-cell adhesion. *Journal of theoretical biology*, 243(1):98-113, 2006.

[3] H. Brezis. Functional analysis, sobolev spaces and partial differential equations, 2011.

[4] A. Buttenschön, T. Hillen, A. Gerisch, and K. J. Painter. A space-jump derivation for non-local models of cell-cell adhesion and non-local chemotaxis. *Journal of mathematical biology*, 76:429-456, 2018.

[5] H. M. Byrne and M. Chaplain. Growth of necrotic tumors in the presence and absence of inhibitors. *Mathematical biosciences*, 135(2):187-216, 1996.

[6] L. G. de Pillis, K. R. Fister, W. Gu, T. Head, K. Maples, T. Neal, A. Murugan, and K. Kozai. Optimal control of mixed immunotherapy and chemotherapy of tumors. *Journal of Biological systems*, 16(01):51-80, 2008.

[7] T. S. Deisboeck and I. D. Couzin. Collective behavior in cancer cell populations. *Bioessays*, 31(2):190-197, 2009.

[8] A. Gerisch and M. A. Chaplain. Mathematical modelling of cancer cell invasion of tissue: local and non-local models and the effect of adhesion. *Journal of theoretical biology*, 250(4):684-704, 2008.

[9] B. Gompertz. Xxiv. on the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. in a letter to francis baily, esq. frs &c. *Philosophical transactions of the Royal Society of London*, (115):513-583, 1825.

[10] A. Keimer and L. Pflug. Existence, uniqueness and regularity results on nonlocal balance laws. *Journal of Differential Equations*, 263(7):4023-4069, 2017.

[11] A. Keimer, L. Pflug, and M. Spinola. Existence, uniqueness and regularity of multi-dimensional nonlocal balance laws with damping. *Journal of Mathematical Analysis and Applications*, 466(1):18-55, 2018.

[12] E. Leschiera, T. Lorenzi, S. Shen, L. Almeida, and C. Audebert. A mathematical model to study the impact of intra-tumour heterogeneity on anti-tumour cd8+ t cell immune response. *Journal of Theoretical Biology*, 538:111028, 2022.

[13] F. Li and B. You. Optimal control for a reaction-diffusion model with tumor-immune interactions. *Communications in Nonlinear Science and Numerical Simulation*, page 108677, 2025.
[14] J.-L. Lions. *Some aspects of the optimal control of distributed parameter systems*. SIAM, 1972.

[15] R. B. Martin. Optimal control drug scheduling of cancer chemotherapy. *Automatica*, 28(6):1113–1123, 1992.

[16] K. J. Painter, J. Bloomfield, J. Sherratt, and A. Gerisch. A nonlocal model for contact attraction and repulsion in heterogeneous cell populations. *Bulletin of mathematical biology*, 77:1132–1165, 2015.

[17] H. Schättler and U. Ledzewicz. Optimal control for mathematical models of cancer therapies. *An application of geometric methods*, 2015.

[18] G. W. Swan. Role of optimal control theory in cancer chemotherapy. *Mathematical biosciences*, 101(2):237–284, 1990.

[19] L. J. Trachette and M. B. Helen. A mathematical model to study the effects of drug resistance and vasculature on the response of solid tumors to chemotherapy. *Journal; Mathematical Biosciences 164 (2000) 17–38*;, (1085), 1999.

[20] F. Tröltzsch. *Optimal control of partial differential equations: theory, methods, and applications*, volume 112. American Mathematical Soc., 2010.

[21] L. Von Bertalanffy. Quantitative laws in metabolism and growth. *The quarterly review of biology*, 32(3):217–231, 1957.

[22] Z. Wu, Z. Luo, and T. Zhang. Analysis and optimal control of a hierarchical age-structured population model in a polluted environment. *Journal of Mathematical Analysis and Applications*, 542(2):128802, 2025.

[23] S. Zheng. *Nonlinear evolution equations*. Chapman and Hall/CRC, 2004.