The diffusion identification in a SIS reaction-diffusion system

1.   Introduction

The mathematical modeling of viral infectious disease transmission is a research area that has received the attention of several researchers in mathematical biology ^{[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]}. Particularly, we refer to ^[16] for a recent summary of challenges in modeling the dynamics of zoonotic infectious diseases. The increasing interest of mathematicians in researching the topic of biology is motivated by several aspects; among them, there is the rapid development of advanced technology and the generation of multidisciplinary teams to understand some complex problems that originated in human body behavior or in the interaction of humans with their environment. The current technology generates numerical datasets that require advanced mathematical tools for their analysis, then there are several mathematical models. The potential contributions of multidisciplinary teams were evidenced in the solution of some specific emergencies generated in the last few years, like the Ebola virus or the COVID-19 pandemic, where the knowledge of mathematicians and epidemiologists was fundamental to the development of health public policy. However, despite the technological development achieved in recent years, there are many aspects that are not observable or measurable by contemporary technology; and, thus, models are necessary to simulate, conjecture; or predict situations where technology is not capable.

Although the current list of works dedicated to the modeling of virus transmission is voluminous, from the mathematical modeling cycle approach; we can distinguish some common steps that are used in the mathematical epidemiology ^[17]: the analysis of the experimental data for the precise disease; the election of the appropriate mathematical framework that models the quantitative data; the mathematical analysis to characterize the behavior of the model; the parameter calibration of the mathematical model; the validation model and improvement of the model if it is necessary, since the process is cyclic. Moreover, in the context of modeling, there are articles that are specialized or are interested in a single step of the cycle; for instance, the articles covering the topic of the well-posedness of the ordinary differential systems that are deduced as mathematical models. In this paper, we are interested in the calibration of a model of the dynamics of two populations during infectious diseases. Then, in order to present the identification problem, we begin by precisely stating the mathematical model.

The epidemiological model considered in this article is based on partial differential equations, where the main variables model the dynamics of two populations: The population of healthy persons who are susceptible to contracting an infection and the population of infected persons who can transmit the disease to the population of healthy ones. More precisely, we consider the following initial boundary value problem:

where $S$ and $I$ are the unknowns of the system, denoting the density of susceptible and infective populations; $\Omega\subset \mathbb{R}^d$ ($d = 1, 2, 3$) is the spatial domain where the total population lives; $\boldsymbol{\nu}$ is the outward normal to the boundary $\partial\Omega$; $\beta$ the transmission rate coefficient and $\gamma$ is the recovery rate coefficient; $d_{ij}, \; i, j = 1, 2$, defined from $\Omega$ to $\mathbb{R}^+$, are the diffusion functions; and $S_0$ and $I_0$ denote the initial densities of susceptible and infective populations. In a broad sense, the system (1.1)–(1.4) is deduced considering the following assumptions: A population living in the spatial domain $\Omega$ is partitioned into two sets of individuals ^[18,19]: Susceptible and infective; the healthy susceptible individuals can contract the disease from cross contacts with infected ones (modeled by the term $\beta SI/(S+I)$); there is no immunity. In that sense, the infected individuals who are recovered can contract the diseases. That fact is modeled by the term $\gamma I$; the spreading of the disease is influenced by the movement of individuals on the domain (modeled by the diffusion terms $d_{ij}$); and we consider that the boundary $\partial\Omega$ is closed to emigration or immigration, which is modeled by the flux boundary condition (1.3).

In order to precisely describe the problem and the main results of the paper, we present the function framework notation and the regularity assumptions on the domain, coefficients, and initial conditions of the problem (1.1)–(1.4). We consider the following function spaces that are standardly used in the analysis of parabolic equations ^[20,21,22]: $C^{k, \alpha}(\overline{\Omega}),$ $k\in\mathbb{N},$ $\alpha \in]0, 1]$ denotes the Hölder $k-$times continuously differentiable functions whose $k^{th}$-partial derivatives are Hölder continuous with exponent $\alpha$; $L^p(\Omega),$ $p\geq 1$; denotes the space of all functions from $\Omega$ to $\mathbb{R}$, which are $p$-integrable in the sense of Lebesgue; $W^{m, p}(\Omega)$ denotes the usual Sobolev spaces of functions that have weak derivatives up to order $m$ and belong to $L^p(\Omega)$; $H^m(\Omega) = W^{m, 2}(\Omega)$ and $C^{\alpha}(\overline{\Omega}) = C^{0, \alpha}(\overline{\Omega})$. Our analysis, in the present study, is conducted by considering the following set of assumptions:

${({\rm{D0}})}$ The spatial environment $\Omega$ is an open bounded and convex set with boundary $\partial\Omega$ of $C^1$ class.

${({\rm{D1}})}$ The functions $S_0$ and $I_0$ defining the initial conditions belonging to $C^{2, \alpha}(\overline{\Omega})$ and satisfying

on $\Omega$, for some positive constant $\phi_0$;

${({\rm{D2}})}$ The coefficients of reaction have the regularity $(\beta, \gamma)\in C^{\alpha}(\overline{\Omega})$ and $(\beta, \gamma)(x) \subseteq[\underline{b}, \overline{b}]\times [\underline{r}, \overline{r}]$ on $\Omega$ for some $\underline{b}, \overline{b}, \underline{r}, \overline{r}\in]0, 1[.$

${({\rm{D3}})}$ The functions $S^{obs}$ and $I^{obs}$ define the observation belonging to $C^{2+\alpha, 1+\alpha/2}(\Omega)$.

We notice that the assumptions (D0)–(D2) are necessary to study the identification problem in the context of classical solutions of the direct problem (1.1)–(1.4), and more weak conditions can be considered to study the problem in the context of weak solutions; see, for instance, ^[23] for the particular case of the identity diffusion matrix.

In this paper, we are interested in the model calibration, or specifically, the main aim is the identification of the diffusion matrix $\mathbb{D} = \mathrm{diag}(d_1, d_2)$ from the observation of susceptible and infective population densities over all spatial domain $\Omega$ in a fixed time $T$, i.e., $S^{obs}, I^{obs}:\Omega\to \mathbb{R}$ are known at time $T$. More precisely, let us consider the functions $\alpha, \gamma, S_0, I_0, S^{obs}$ and $I^{obs}$ from $\Omega$ to $\mathbb{R}^+$ are known, and we want to determine $\mathbb{D}$ the solution of the following constrained optimization problem

where $J$ and $\mathcal{U}_{ad}(\Omega)$ are defined as follows

and $\Gamma > 0$ is an appropriate constant. Hereinafter, we consider the notation $\|\mathbb{D}\|^2_{L^2(\Omega)^2} = \|d_1\|^2_{L^2(\Omega)}+\|d_2\|^2_{L^2(\Omega)}$. We note that the functional $J$ defined on (1.7) is the more pertinent for the determination of $\mathbb{D}$, since the first and second terms of $J$ are a comparison of the state solution profiles $(S(\cdot, T), I(\cdot, T))$ and the observation $(S^{obs}, I^{obs})$ in the $L^2$ norm and the third term is a regularization term, where $\Gamma$ should be appropriately selected in order to get a unique solution of the optimization problem.

SIR (susceptible–infectious–removed) SIS (susceptible–infectious–susceptible)

The analysis for the calibration of compartmental models (susceptible–infectious–susceptible and susceptible–infectious–removed shortly as SIS and SIR, respectively) was recently developed ^{[17,19,23,24,25,26,27]}. We observe that in all of those works, the authors identify the reaction coefficients; for instance, in ^[19] the authors get results for the identification of reaction term coefficients in the one-dimensional spatial domain ($d = 1$), and it is extended to the higher dimensions ($d\ge 2$) in ^[17,24]. Moreover, in those works, the matrix modeling the diffusion is the identity matrix. However, to the best of our knowledge, the identification of matrix diffusion in epidemiological compartmental models has not yet been conducted. However, we must recognize that there are some works on the identification of the diffusion matrix in linear elliptic and parabolic problems whose results are not directly extensible to matrix diffusion identification in nonlinear systems of reaction-diffusion ^{[28,29,30,31,32]}.

The main results, which are the contributions of this paper, are given by the following five results: (ⅰ) The introduction of the necessary conditions to establish the existence and uniqueness of a positive solution to the direct problem (1.1)–(1.4) (see Section 2); (ⅱ) the existence of optimal solutions for (1.5)–(1.9) (see Section 3); (ⅲ) the introduction of an adjoint system with classical bounded solution (see Section 4); (ⅳ) the definition of a first-order optimality condition that characterizes the optimal solution in terms of direct and adjoint state solutions (see Section 5); and (ⅴ) a local uniqueness of identification problem (see Section 6). Furthermore, we present two numerical examples on and state some main conclusions (see Sections 7 and 8, respectively).

2.   Well-posedness of the state system (1.1)–(1.4)

Theorem 2.1. Consider that the hypotheses (D0)–(D2) are satisfied. If $(d_1, d_2)\in C^{\alpha}(\overline{\Omega})^2$, there is a unique positive pair of functions $(S, I)\in C^{2+\alpha, 1+\alpha/2}(\overline{Q}_T)^2$ that satisfies the direct problem defined by the initial boundary value problem (1.1)–(1.4), which admits a unique positive classical solution $(S, I)$. Moreover, $S$ and $I$ are bounded on $\overline{Q}_T$, i.e., the estimate

is satisfied for any given $T\in \mathbb{R}^+$.

Proof. If we assume the existence of the solution of (1.1)–(1.4), we deduce some a priori estimates. We can prove the nonnegativity behavior of $S$ and $T$ by applying the maximum principle. From (1.1), (1.2), the positivity of $S$ and $T$, the relation $S/(S+I) < 1$; and the bounded behavior of reaction coefficients, we deduce that

Moreover, we can observe that $(S_{\max}, I_{\max})$ is a supersolution of the following linear system

Thus, the upper a priori estimate

can be deduced by applying the well-known comparison principle for parabolic equations.

We can follow the local existence of classical solutions of (1.1)–(1.4) by the standard results given in ^[33,34,35] and we can deduce the Hölder regularity of the local solution by modifying appropriately the arguments used in ^[18]. Thus, we have that there is a pair of nonnegativity functions $(\hat{S}, \hat{I})(x, t)$ that are the local solutions of (1.1)–(1.4); or, equivalently, there is $T_{\max} > 0$ (the maximal existence time), such that $(\hat{S}, \hat{I})(x, t)$ is the following initial boundary value problem:

Consequently, the proof of existence and uniqueness of the global solution is reduced to guaranteeing the L$^\infty$ estimations of $\hat{S}$ and $\hat{T}$ and applying similar arguments to those given in ^[18,36] (see also Theorem 2 in ^[37] and Lemma 1.1 in ^[38]).

We observe that the positive constant $p_0$ defined on the relation (1.6) of ^[37] should be selected such that $p_0 > d\max\{0, 1\}/2 > 3/2$ [37, Theorem 1], then to obtain the L$^\infty$ estimations is enough to obtain estimates in L$^{p_0}$ for some $p_0 > 3/2$. Indeed, we select $p_0 = 2$; and we derive L$_2$ estimations of $(\hat{S}, \hat{T})$. Multiplying (2.3) by $\hat{S}$, integrating on $\Omega$; and (2.2), we have that

Similarly, multiplying (2.4) by $\hat{I}$; integrating on $\Omega$, and using the fact that $\hat{S}/(\hat{S}+\hat{I})\le 1$, we have that

The estimates (2.7) and (2.8) implies that

and, consequently, with the application of [37, Theorem 1], we deduce the existence and uniqueness of the global solution and, particularly, the estimate (2.1) is satisfied. □

3.   Existence of solutions for the constrained optimization problem (1.5)–(1.9)

Theorem 3.1. Consider the assumptions (D0)–(D3) are satisfied, then the optimization problem (1.5)–(1.9) has at least one solution.

Proof. We note that the admissible set is not empty and $\mathcal{J}(\mathbb{D})$ is bounded for any $\mathbb{D}\in \mathcal{U}_{ad}(\Omega).$ The first assertion, i.e., $\mathcal{U}_{ad}(\Omega)\not = \emptyset$, follows by considering the diffusion matrix $\mathbb{D}(x) = \mathrm{diag}(\underline{\delta_1}+\overline{\delta_1}, \underline{\delta_2}+\overline{\delta_2})/2\in\mathcal{U}_{ad}(\Omega)$. Meanwhile, we can prove that the cost function $\mathcal{J}$ is bounded by analyzing the boundedness of each term: The first two terms are bounded as consequence of the bounded behavior of the direct problem as result of Theorem 2.1, and the regularity of the observation functions is given on hypothesis (D3); the third term is bounded as consequence of the fact that $\mathbb{D}\in \mathcal{U}_{ad}(\Omega)$ and the definition of the admissible set. Consequently, we can deduce the existence of as a minimizing sequence of $\mathcal{J}$, where $\mathcal{M}$ is a bounded and closed set of $H^{{|[d/2]|} +1}(\Omega)^2.$ We observe that the following compact embedding $H^{{|[d/2]|} +1}(\Omega)\subset C^\alpha (\Omega)$ is satisfied for all $\alpha\in]0, 1/2]$ and the convexity of $\Omega$ is assumed on (D0). This kind of inclusion is the consequence of two results: $H^{{|[d/2]|} +1}(\Omega)$ is continuous embedding in $C^{1/2}(\Omega)$ (see Theorem 6 [39, pp 270]), and $C^{1/2}(\Omega)$ is compact embedding in $C^{\alpha} (\Omega)$ for all $\alpha\in]0, 1/2]$, and $\Omega$ is a convex set (see Theorem 1.3.1 [40, pp 11]). Thus, clearly, $H^{{|[d/2]|} +1}(\Omega)\subset C^{1/2}(\Omega)\subset C^{\alpha} (\Omega)$ for all $\alpha\in]0, 1/2]$ implies that the embedding $H^{{|[d/2]|} +1}(\Omega)$ in $C^\alpha (\Omega)$ is compact for all $\alpha\in]0, 1/2]$, and $\Omega$ is a convex set.

The compact embedding of $H^{{|[d/2]|} +1}(\Omega)$ in $C^\alpha (\Omega)$ for $\alpha\in]0, 1/2]$ and $\Omega$ convex, implies that $\{\mathbb{D}_n\}$ is bounded in the strong topology of $C^{\alpha}(\overline{\Omega})^2$ for all $\alpha\in]0, 1/2],$ since

where $C$ is independent of $d_1, d_2$ and $n$. Here, we remark that the righthand side is bounded by the fact that with $\mathcal{M}$ as a bounded and closed set of $H^{{|[d/2]|} +1}(\Omega)^2.$

Let us consider the notation $(S_n, I_n)$ to the solution of the direct problem (1.1)–(1.4) corresponding to $\mathbb{D}_n$, then, by considering the fact that $\{\mathbb{D}_n\}\in C^{\alpha}(\overline{\Omega})^2$ for all $\alpha\in]0, 1/2]$, by Theorem 2.1, we have that $(S_n, I_n)\in C^{2+\alpha, 1+\frac{\alpha}{2}}(\overline{Q}_T)^2$. Also $\{(S_n, I_n)\}$ is a bounded sequence in the strong topology of $C^{2+\alpha, 1+\frac{\alpha}{2}}(\overline{Q}_T)^2$ for all $\alpha\in]0, 1/2]$.

The boundedness of the sequence $\{(\mathbb{D}_n, S_n, I_n)\}$, implies that there exists $(\overline{\mathbb{D}}, \overline{S}, \overline{T})$ such that

and uniformly convergent subsequences, which are again labeled by $\{\mathbb{D}_n\}$ and $\{(S_n, I_n)\}$; to be precise

Moreover, it is straightforward to deduce that $(\overline{S}, \overline{I})$ is the solution of (1.1)–(1.4) when the diffusion matrix is given by $\overline{\mathbb{D}}.$ Hence, using the definition of the minimizing sequence, the weak lower-semicontinuity of the $L^2$ norm, and the Lebesgue's dominated convergence theorem, we get that

Thus, $\overline{\mathbb{D}}$ is a solution of (1.5)–(1.9). □

4.   Adjoint system

In order to deduce the adjoint system, we adapt the formal calculus of the adjoint equation for scalar strongly parabolic equations in ^[41,42]. Let us consider $\mathcal{L}$, the Lagrangian associated to the optimization problem (1.5)–(1.9), defined as follows

where $E_1$ and $E_2$ are the weak formulations of (1.1) and (1.2), respectively. More precisely

for $p$ and $q$ the test functions.

Let $\overline{\mathbb{D}}$ be a solution of optimization problem (1.5)–(1.9) and $(\overline{S}, \overline{I})$ be the solution of the forward problem (1.1)–(1.4) with $\overline{\mathbb{D}}$ instead of $\mathbb{D}$. By a formal calculus of the derivative of $\mathcal{L}$ with respect to $d_1$ and $d_2$ and introducing the test functions $(p, q)$, such that the derivatives of the state variables $(\overline{S}, \overline{I})$ with respect to $d_1$ and $d_2$ are vanished, we get that the functions $(p, q)$ are obtained as the solution of the backward boundary value problem:

The system (4.2)–(4.5) is called the adjoint system to (1.1)–(1.4).

Theorem 4.1. Consider that $\Omega, S_0, I_0, S^{obs}, I^{obs}, \beta$ and $\gamma$; satisfy the assumptions of Theorem 3.1. Moreover, consider that $\overline{\mathbb{D}}\in U_{ad}$ is a solution of (1.5)–(1.9); and $(\overline{S}, \overline{I})$ is the solution of the direct problem (1.1)–(1.4) with $\overline{\mathbb{D}}$ instead of $\mathbb{D}$, then, the solution of (4.2)–(4.5) satisfies the following estimates

for $t\in [0, T]$. Here, $C$ denotes some positive generic constant.

Proof. If we introduce the change of the time variable by the following relation $\tau = T-t$ for $t\in [0, T]$ and the unknowns of the direct problem and the adjoint system by the identity $(w_1, w_2, S^{*}, I^{*})(x, \tau) = (p_1, p_2, \overline{S}, \overline{I})(x, T-\tau)$, we can rewrite the adjoint system (4.2)–(4.5) as the following initial boundary value problem

Next, by applying the standard arguments of energy and regularity of solutions for linear parabolic equations, we get the desired estimates (4.6)–(4.9), see ^[19] for $d = 1$ and ^[24] for $d\ge 1$ for the case of the identity matrix diffusion. □

5.   Necessary optimality conditions.

Theorem 5.1. Consider that $\overline{\mathbb{D}}$ is a solution of the optimization problem (1.5)–(1.9), $(\overline{S}, \overline{I})$ is the solution of the direct problem (1.1)–(1.4) with $\overline{\mathbb{D}}$ instead of $\mathbb{D},$ and $(p, q)$ is the solution of the adjoint system (4.2)–(4.5), then, the inequality

is satisfied and defines the first-order optimality condition.

Proof. Let us consider the arbitrary diffusion $\hat{\mathbb{D}}\in U_{ad}(\Omega)$, then we define the notation

where $(S^{\varepsilon}, I^{\varepsilon})$ is the solution of (1.1)–(1.4) with $\mathbb{D}^{\varepsilon}$ instead of $\mathbb{D}$. The fact that that $\overline{\mathbb{D}}$ is an optimal solution of (1.5)–(1.9), by taking the Fréchet derivative of $J_{\varepsilon}$, we deduce that the following inequality

is satisfied. Here, $\partial_{\varepsilon}S^{\varepsilon}$ and $\partial_{\varepsilon}I^{\varepsilon}$ for $\varepsilon = 0$ are the sensitivities of solutions for (1.1)–(1.4), with respect to the $\varepsilon$-perturbations of $\overline{\mathbb{D}}$.

The calculus of the sensitivities $(\partial_{\varepsilon}S^{\varepsilon}, \partial_{\varepsilon}I^{\varepsilon})$ when $\varepsilon\to 0$ is developed by considering the SIS systems of the form (1.1)–(1.4), $(S^{\varepsilon}, I^{\varepsilon})$ and $(\overline{S}, \overline{I})$, then letting $\epsilon\to 0$. More precisely, we have that $(S^{\varepsilon}, I^{\varepsilon})$ and $(\overline{S}, \overline{I})$ are solutions of the following initial boundary value problems

and

respectively. Subtracting the system (5.7)–(5.10) from (5.3)–(5.6), dividing by $\varepsilon$ and using the notation $\left(z_{1}^{\varepsilon}, z_{2}^{\varepsilon} \right) = \varepsilon^{-1}\left(S^{\varepsilon}-\overline{S}, I^{\varepsilon}-\overline{I} \right)$, we deduce the initial boundary value problem

Let us consider that $(z_{1}, z_{2})$ is the limit of $(z_{1}^{\varepsilon}, z_{2}^{\varepsilon})$ when $\varepsilon\to 0$, from (5.11)–(5.14); we deduce straightforward answer that $(z_{1}, z_{2})$ is a solution of the following system

We remark that, in the context of optimization with partial differential equation constraints, the system (5.15)–(5.18) is called the sensitivity system for (1.1)–(1.4).

Using the sensitivity system (5.15)–(5.18), we observe that the relation (5.2) can be rewritten as follows

Moreover, we notice two facts: First

and second, by easy algebraic computations, from the systems (4.2)–(4.5) and (5.15)–(5.18), we can deduce the following identity

which implies that

by integration on $Q_{T}$. Thus, the relations (5.21) and (5.20) implies that

Hence, we can conclude the proof of (5.1) by replacing (5.22) in the first term of (5.19). □

6.   A result of local uniqueness

Theorem 6.1. Let us consider that (D0)–(D3) is satisfied and let us consider the quotient set $\mathcal{U}(\Omega)/\sim$ where the equivalence relation $\sim$ is defined as follows

then, there exists $\Gamma^*\in\mathbb{R}^+$ such that the solution of the optimization problem (1.5)–(1.9) is uniquely defined (up an additive constant) on the quotient set $\mathcal{U}(\Omega)/\sim$ for any regularization parameter $\Gamma > \Gamma^*$.

Proof. Let us consider that $\mathbb{D}, \hat{\mathbb{D}}\in \mathcal{U}(\Omega)/\sim$ are two solutions of (1.5)–(1.9). Moreover, let us consider that the sets of functions $\{S, I, p, q\}$ and $\{\hat{S}, \hat{I}, \hat{p}, \hat{q}\}$ are solutions to the systems (1.1)–(1.4) and (4.2)–(4.5) with diffusion matrices $\mathbb{D}$ and $\hat{\mathbb{D}},$ respectively. From Theorem 5.1 and the hypothesis that $\mathbb{D}$ and $\hat{\mathbb{D}}$ are solutions of (1.5)–(1.9), we have that the following inequalities

are satisfied. If we choose the particular cases $\overline{\overline{{\mathbb{D}}}} = \hat{\mathbb{D}}$ in (6.2) and $\underline{\underline{{\mathbb{D}}}} = \mathbb{D}$ in (6.3), and then add both inequalities, we get

We observe that, by integrating by parts two times, applying the Theorems (2.1) and (4.1), and using the fact that $\mathbb{D}, \hat{\mathbb{D}}\in \mathcal{U}(\Omega)/\sim$, we can bound the righthand side of (6.4), as follows

with

Thus, from (6.5) and (6.4), we deduce the desired uniqueness result. □

7.   Numerical examples

In this section, we consider two numerical examples for the one-dimensional case where the identification is developed from observations that are constructed by considering synthetic data as observation of state variables. We begin by stating precisely that we introduce a small modification of the notation introduced previously. The diffusion coefficients $d_1$ and $d_2$ depend on a finite number of parameters denoted by ${\bf e} = (e_1, \ldots, e_k)\in \mathbb{R}^k$, which is explicitly denoted by $d_i(x) = d_i(x; {\bf e})$ for $i = 1, 2.$ The system (1.1)–(1.4) is modified by considering the mass action $\beta(x) SI$ instead of $\beta(x) SI/(S+I)$. We consider that kind of modification in the direct problem, since our aim is to apply the unconditionally stable Implicit-Explicit (IMEX) numerical method introduced by ^[43]. To be precise, the direct problem considered for the numerical simulations is given by the following initial boundary value problem:

where $\Omega = ]0, 1[, \partial\Omega = \{0, 1\}$ and $\Gamma_T = \{0, 1\}\times [0, T]$. Concerning to the discretization of $Q_T$, we select $M, N\in \mathbb{N}$ such that the discretization of $\Omega$ is given by $x_k = k\Delta x$ for $k = 0, \ldots, M+1$ with $\Delta x = 1/(M+1)$, and the discretization of $[0, T]$ is given by $t_n = n\Delta t$ for $n = 0, \ldots, N$ with $\Delta t = 1/N$. The approximation of a given function $H:\Omega\times \mathbb{R}_+\to \mathbb{R}$ at $(x_k, t_n)$ is denoted by $H^n_k$. Adapting the numerical discretization throughout a finite differences scheme introduced in ^[43], we deduce that the approximation of the initial boundary value problem (7.1)–(7.5) is given by

We remark that the numerical method considered in ^[43] is developed and analized when $d_1$ and $d_2$ are constants. However, by straightforward adaptation of the arguments given in ^[43], we can deduce that the implicit–explicit numerical scheme (7.6)–(7.9) has several properties, basically preserves the biological meaning (such as positivity), and is unconditionally convergent.

For discretization of the inverse problem (1.5) and (1.6), we begin by considering the discretized cost function

We observe that in (7.10) we have omitted the regularization term, i.e. $\Gamma = 0,$ then, in our numerical example, we consider that the inverse problem (1.5) and (1.6), is replaced by the following parameter identification problem

In both numerical examples, we solve the optimization problem using the optimset routine of Matlab.

7.1.   Example 1: Identification of a constant diffusion

We select the following coefficients on the reaction term $\beta(x) = 0.000284535$ and $\gamma(x) = 0.144$ and the initial condition $(S_0, I_0)(x) = (x, 2-x)/2$. We consider that ${\bf e} = (e_1, e_2)$, $d_1(x; {\bf e}) = e_1$ and $d_2(x; {\bf e}) = e_2$. We construct the observation profile at $T = 0.6$ by considering a numerical simulation of the direct problem with ${\bf e}^{obs} = (0.5, 0.5)$, $M = 200$ and $N = 100000$ (i.e., $\Delta x = 5E-3$ and $\Delta x = 6E-6$). The state simulation on $Q_T$ is shown on Figure 1(a), (b). We consider the initial guess ${\bf e}^{obs} = (0.1, 0.1)$ and get that the identified parameters are ${\bf e}^{\infty} = (0.52294, 0.55149)$. The numerical identification is developed by considering $M = 100$ and $N = 1000$ or, equivalently, $\Delta x = 1.0E-2$ and $\Delta x = 5.\overline{9}E-4$. The comparison of the observed, identified and initial guess profiles are shown in Figure 1(c)–(f).

7.2.   Example 2: identification of a quadratic diffusion function

In this numerical example, we consider $\beta(x) = 0.000284535$ and $\gamma(x) = 0.144$ and the initial condition $S_0(x) = 5$ and

The parameters to identify are given by ${\bf e}\in \mathbb{R}^4$, such that the diffusion functions are of the following parametric form

The independent terms in $d_1$ and $d_2$ are fixed to prevent the degeneration of the diffusion function. The observation profiles $T = 0.6$ are constructed by a numerical simulation of the direct problem with ${\bf e}^{obs} = (0.5, 0.5, 0.5, 0.5)$, $M = 200$ ($\Delta x = 5E-3$) and $N = 100000$ ($\Delta x = 6E-6$), which are shown in Figure 2(a), (b). The initial guess and identified parameters are given by ${\bf e}^{obs} = (0.1, 0.1, 0.1, 0.1)$ and ${\bf e}^{\infty} = (0.75143, 0.42256, 0.95842, 0.21146)$, respectively. For the identification, we assume that the discretization $M = 100$ ($\Delta x = 1.0E-2$) and $N = 1000$ ($\Delta x = 5.\overline{9}E-4$). The comparison of observed, identified and initial guess profiles and diffusion functions are shown in Figure 2(c)–(f).

From Figure 2(c), (e), we observe that the profile $S_{{\bf e}^{\infty}}(\cdot, 0.6)$ fits the observation data. However, $d_1(x, {\bf e}^{\infty})$ is close to $d_1(x, {\bf e}^{obs})$, but we conjecture that it can be improved by incorporating the regularization term on $J_{\Delta}$. A similar behavior is observed from Figure 2(d), (f) for the cases of $I_{{\bf e}^{\infty}}(\cdot, 0.6)$ and $d_2(x, {\bf e}^{\infty})$.

8.   Conclusions

In this paper, we have introduced the functional framework to develop the identification of the diffusion matrix in a reaction-diffusion system arising from the modeling of the spread dynamics of virus propagation. We considered that the disease occurs in a spatially distributed population between two classes of individuals: The susceptible class, formed by the individuals who can catch the disease; and the infective class, formed by the individuals who are infected and can transmit the disease. The reaction-diffusion model was deduced by assuming that there are no vital dynamics, there is no migration or immigration during the epidemic disease propagation, and the coefficients (diffusion, transmission rate, and recovery rate) of the model are functions that depend on the spatial position. The diffusion matrix identification was developed by assuming that the susceptible and infective populations are known at a fixed time. Thus, we have formulated the inverse problem as an optimal control problem, where the cost function to minimize is the least squares cost function and a regularization term, and the optimization constraints are the SIS reaction-diffusion model.

We have proved that the mathematical model is well-posed and has a global positive solution in the context of strong solutions in Hölder spaces when the initial conditions and coefficients are of Hölder class. We have demonstrated that there exists at least one solution to the identification problem, and the solution is unique under the assumption that the regularization parameter is large enough. Furthermore, we remarked that the uniqueness of the optimal control problem is deduced from an appropriate quotient set of the admissible set. Moreover, we have introduced a necessary optimal for the optimal control problem.

We observed that our results define the appropriate framework to develop the numerical identification from available experimental data for a concrete epidemic propagation and even numerical analysis like convergence. Here, we can remark that in a practical epidemiology phenomenon, the observation data of profiles in all spatial domain are not usual, and the typical situation is a profile at a fixed point in the domain and during a time interval. The cost function considered in this paper must be modified to develop identification from experimental data. Moreover, in the case of parameter identification from laboratory or epidemic data, we should develop a study of the noise, and we can consider the recent uncertainty concepts ^[44].

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

A. C. and F. H. thank the partial support of Universidad del Bío-Bío (Chile) through the research projects: 2120436 IF/R, research project INES I+D 22-14, and FAPEI, and the research project of the Postdoctoral Program as a part of the project "Instalación del Plan Plurianual UBB 2016-2020"; National Agency for Research and Development, ANID-Chile, through FONDECYT project 1230560; Universidad Tecnológica Metropolitana through the project supported by the Competition for Research Regular Projects, year 2020, Code LPR20-06, Universidad Tecnológica Metropolitana. A. T. and I. H. acknowledge the partial support of "ayudantes de investigación" at Universidad del Bío-Bío (Chile) and ANID (Chile) through the program "Becas de Doctorado".

Conflict of interest

The authors declare there is no conflict of interest.