Mathematical analysis, forecasting and optimal control of HIV/AIDS spatiotemporal transmission with a reaction diffusion SICA model

1.   Introduction

The human immunodeficiency virus (HIV) causes millions of deaths to humans worldwide, being one of the most infectious and deadly virus ^[10]. The deterministic SICA model was introduced by Silva and Torres in 2015, as a sub-model of a general Tuberculosis and HIV/AIDS (acquired immunodeficiency syndrome) co-infection problem ^[11]. After that, it has been extensively used to investigate HIV/AIDS, in different settings and contexts, using fractional-order derivatives ^[14], stochasticity ^[1] and discrete-time operators ^[17], and adjusted to different HIV/AIDS epidemics, as those in Cape Verde ^[12] and Morocco ^[8].

One of the fundamental characteristics of SICA modeling is that it provides adequate but simple mathematical models that help to characterize and understand some of the essential epidemiological factors leading to the spreed of the AIDS disease. In such models, the susceptible population $S$ is nourished by the recruitment of individuals into the population at a rate $\lambda$. All individuals are exposed to natural death, at a constant rate $\mu$. Individuals $S$ are susceptible to HIV infection from an effective contact with an individual carrying the HIV, at the rate $\dfrac{\beta}{N}\left(I+\eta_C C+\eta_A A\right)$, where $I$, $C$ and $A$ denote, respectively, the infected, chronic (under treatment) and AIDS individuals, $N$ represents the total number of individuals in the population under study, that is, $N$ is the sum of $S$, $I$, $C$ and $A$ individuals, and $\beta$, $\eta_C$ and $\eta_A$ are parameters that depend on the particular situation under study. For a survey on SICA models for HIV transmission, showing that they provide a good framework for interventions and strategies to fight against the transmission of the HIV/AIDS epidemic, we refer the reader to ^[15].

It is well known that reaction-diffusion equations are commonly used to model a variety of physical and biological phenomena ^{[2,4,6,16,19,21]}. Such equations describe how the concentration or density distributed in space varies under the influence of two processes: (i) local interactions of species and (ii) diffusion, which causes the spread of species in space. Recently, reaction-diffusion equations have been used by many authors in epidemiology as well as virology, see, e.g., ^[20], where a mathematical model is proposed to simulate the hepatitis B virus infection with spatial dependence, or the non-theoretical reviews ^[3,5]: in ^[3], host-pathogen interactions are described by different temporal and spatial scales, while ^[5] covers bioinformatics workflows and tools for the routine detection of the SARS-CoV-2 infection. Here we propose, for the first time in the literature, to use SICA modeling with $S$, $I$, $C$ and $A$ (thus, also $N$) as functions of both time $t$ and space $x$. The spatial effect plays a crucial role in the spread of the virus. In order to well describe this phenomenon, we incorporate terms that model the spatial diffusion in each compartment, by adding $\Delta S$, $\Delta I$, $\Delta C$ and $\Delta A$ in the classical SICA model system. By taking into account the spatiotemporal diffusion allow us not to neglect a good part of compartments' inputs-outputs.

The paper is organized as follows. We begin with some preliminaries on the physical interpretation of the Laplacian in Section 2. The spatiotemporal SICA model is then introduced in Section 3 and its mathematical analysis is given in Section 4 where, by using semigroup theory ^[9,18], we prove existence and uniqueness of a strong nonnegative solution to the system (see Theorem 4.1). In Section 5, we show some numerical examples that motivate us to consider optimal control. An optimal control problem is then formulated and existence of a solution is established (see Theorem 5.1). Next, we obtain in Section 6 a set of necessary optimality conditions that characterize the optimal solution. We end with Section 7 of conclusions, pointing also some future directions of research.

2.   Preliminaries: interpretation of the Laplacian

Let $\nabla^2$ be the Laplacian in two dimensions expressed by

Suppose that, at a point $O$, taken as the origin of the system of axises $Oxy$, a field $f$ takes the value $f_0$. Consider an elementary square with side $a$ whose edges are parallel to the coordinate axises and whose center merges with the origin $O$. The average value of $f$ in this elementary cube, that is, the mean value of $f$ in the neighborhood of the point $O$, is given by the expression

where the two integrations relate to the rectangle $C = [-\frac{a}{2}, \frac {a}{2}]^2$. At an arbitrary point $P(x, y)$ in the neighborhood of $O = (0, 0)$, we develop $f$ in Taylor–Maclaurin series. Thus,

On one hand, the odd functions in this expression provide, by integration from $-\frac{a}{2}$ to $\frac {a}{2}$, a zero contribution to $\overline{f}$. For example,

On the other hand, each even function provide a contribution of $\frac{a^4}{12}$. For example,

Using the Fubini–Tonnelli theorem, we get

We deduce that

and

As the point $O$ has been chosen arbitrarily, we can assimilate it to the current point $P$ and drop the index $0$. Therefore, we obtain the expression

the interpretation of which is immediate: the quantity $\nabla^2 f$ is approximately proportional to the difference $\overline {f} - f$. The constant of proportionality is worth $\frac{24}{a^4}$ in Cartesian axises. In other words, the quantity $\nabla^2 f$ is a measure of the difference between the value of $f$ at any point $P$ and the mean value $\overline{f}$ in the neighborhood of point $P$.

3.   The spatiotemporal mathematical SICA model

In ^[12], Silva and Torres proposed the following epidemic SICA model:

The limitation of the temporal dynamical system (3.1) to give a good description of the spread of the virus in the space is obvious. To bridge this gap, we suggest to use of the Laplacian operator as interpreted in Section 2. In concrete, we extend the deterministic epidemic SICA model (3.1) as follows:

where $\Delta$ is the Laplacian in the two-dimensional space $(t, x)$ and $u:[0;T]\times \Omega\longrightarrow [0;1[$ is a control that permits to diminish the number of infected individuals and to increase that of susceptible by devoting some special treatment to the most affected persons. The description of the parameters of model (3.2) is summarized in Table 1.

4.   Existence and uniqueness of a strong nonnegative solution

In order to prove existence and uniqueness of a strong solution to system (3.2), we define some tools. Consider the Hilbert spaces $H(\Omega) = (L_2(\Omega))^4$, $H^1(\Omega) = \left\{u\in L_2(\Omega):\; \dfrac{\partial u}{\partial x}\in L_2(\Omega)\; {\rm{and}}\; \dfrac{\partial u}{\partial y}\in L_2(\Omega)\right\}$ and $H^2(\Omega) = \left\{u\in H^1(\Omega):\; \dfrac{\partial^2 u}{\partial x^2}, \dfrac{\partial^2 u}{\partial y^2}, \dfrac{\partial^2 u}{\partial x\partial y}, \dfrac{\partial^2 u}{\partial y\partial x} \in L_2(\Omega)\right\}$. Let $L^2(0, T;H^2(\Omega))$ be the space of all strongly measurable functions $v:[0, T]\longmapsto H^2(\Omega)$ such that

and $L^{\infty}(0, T; H^1(\Omega))$ be the set of all functions $v:[0, T]\longmapsto H^1(\Omega)$ verifying

The norm in $L^{\infty}(0, T; H^1(\Omega))$ is defined by

Our model is equivalent to

where $z = (z_1, z_2, z_3, z_4) = (S, I, C, A)$ and $g = (g_1, g_2, g_3, g_4)$ is defined by

For all $i\in \{1, 2, 3, 4\}$,

Let $A$ denote the linear operator defined from $D(A)\subset H(\Omega)$ to $H(\Omega)$ by

with

and $U_{ad}$ be the admissible control set defined by

with $Q = [0, T]\times \Omega$ and $\Omega$ a bounded domain in $\mathbb{R}^2$ with smooth boundary $\partial \Omega$.

To obtain our next result, we employ semi-group theory ^[18] to prove existence and uniqueness of a global nonnegative solution to the considered system.

Theorem 1. Let $\Omega$ be a bounded domain from $\mathbb{R}^2$ with a boundary of class $C^{2+\alpha}$, $\alpha > 0$. For nonnegative parameters of the spatiotemporal SICA model (3.2), $u\in U_{ad}$, $z^0\in D(A)$ and $z^0_i\geq 0$ on $\Omega$, $i = 1, 2, 3, 4$, the system (3.2) has a unique (global) strong nonnegative solution $z \in W^{1, 2}([0, T];H(\Omega))$ such that

Additionally, there exists $C > 0$, independent of $u$ and of the corresponding solution $z$, such that for all $t \in [0, T]$ and all $i\in \{1, 2, 3, 4\}$ one has

Proof. Because the Laplacian operator $\Delta$ is dissipating, self-adjoint, and generates a $C_0-$ semigroup of contractions on $H(\Omega)$, it is clear that function $g = (g_1, g_2, g_3, g_4)$ becomes Lipschitz continuous in $z = (z_1, z_2, z_3, z_4)$ uniformly with respect to $t \in [0, T]$. Therefore, the problem admits a unique strong solution $z$. Let us now show that for all $i\in \{1, 2, 3, 4\}$, $z_i\in L_{\infty}(Q)$. Indeed, set $k = \max\left\{\|g_i\|_{L_{\infty}(Q}), \|z_i^0\|_{L_{\infty}(\Omega)}\; :\; i\in \{1, 2, 3, 4\} \right\}$ and let

Then,

Let $i\in \{1, 2, 3, 4\}$. There exists an infinitesimal semigroup $\Gamma(t)$ associated to the operator $d_i\Delta$ such that

We deduce that $U_i (t, x)\leq 0$ and so $z_i\leq kt+\|z_i^0\|_{L^{\infty}(\Omega)}$.

Consider $V_i (t, x) = z_i (t, x)+kt+\|z_i^0\|_{L^{\infty}(\Omega)}$. Upon differentiation, we get

The strong solution of the above equation is

Then, $V_i (t, x)\geq 0$ and so $z_i\geq -kt-\|z_i^0\|_{L^{\infty}(\Omega)}$. Consequently, $|z_i(t, x, )|\leq kt+\|z_i^0\|_{L^{\infty}(\Omega)}$, which implies that $z_i\in L_{\infty}(Q)$.

Now, we proceed by proving that $z_i\in L_{\infty}\left(0, T;H^1(\Omega)\right)$ for all $i\in \{1, 2, 3, 4\}$. Indeed, let $i\in \{1, 2, 3, 4\}$. From equality

we obtain that

From Green's formula, we get

Since $g_i\in L^2(Q)$, $z^0_i\in L^2(Q)$ and $z_i, z_i^0\in L_{\infty}(Q)$, we obtain that $z_i\in L_{\infty}\left(0;T; H^1(\Omega))\right)$.

Finally, using the same arguments as for the Field–Noyes equations in [16,Example 4], we deduce that the solution $(z_1, z_2, z_3, z_4)$ is nonnegative. Consider the set

and the convex functions $G_i$ defined on $\Sigma$ by $G_i(z_1, z_2, z_3, z_4) = -z_i$. One can see that

According to [16,Theorem 14.14], the region $\Sigma$ is positively invariant and the result follows.

5.   Existence of an optimal control

To motivate the interest on optimal control, we begin by showing some numerical simulations of our spatiotemporal SICA model (3.2). For details on the simulation method, tool and used code, see Appendix A.

We have considered the values for the parameters as given in Table 2, which were borrowed from ^[12].

Then, the dynamics without control, that is, with $u \equiv 0$ in (3.2), is given in Figure 1.

In contrast, dynamics in the presence of a control are given in Figures 2 and 3. We conclude that the evolution of the system related with the absence of control differs totally to those in presence of controls. Indeed, Figure 1 shows that in absence of the control the density of the infected individuals increases while in the presence of a control (Figures 2 and 3) it clearly decreases. The question of how to choose the control along time, in an optimal way, is therefore a natural one.

Motivated by ^[13], our aim is to minimize the sum of the density of infected individuals and the cost of the treatment program. Mathematically, the problem we consider here is to minimize the objective functional

subject to the control system (3.2) and where the admissible control set $U_{ad}$ is defined as in (4.2).

Theorem 2. Under the conditions of Theorem 1, our optimalcontrol problem admits a solution $(z^*, u^*)$.

Proof. The proof is divided into three steps.

Step 1: Existence of a minimizing sequence $(z^n, u_n)$. The infimum of the objective function on the set of admissible controls is ensured by the positivity of $J$. Assume that $J^* = \inf_{u\in U_{ad}}J(z, u)$. Let $\{u_n\}\subset U_{ad}$ be a minimizing sequence such that $\lim\limits_{n\rightarrow +\infty}J(z^n, u_n) = J^*$, where $(z^n_1, z^n_2, z^n_3, z^n_4)$ is the solution of the system corresponding to the control $u_n$. Subsequently,

where $\dfrac{\partial z^n_1}{\partial \eta} = \dfrac{\partial z^n_2}{\partial \eta} = \dfrac{\partial z^n_3}{\partial \eta} = \dfrac{\partial z^n_4}{\partial \eta} = 0$ on $Q$.

Step 2: Convergence of the minimizing sequence $(z^n, u_n)$ to $(z^*, u^*)$. Let $i\in\{1, 2, 3, 4\}$. Note that $z_i^n (t, x)$ is compact in $L^2(\Omega)$ from the fact that $H^1(\Omega)$ is compactly embedded in $L^2(\Omega)$. In order to apply the Ascoli–Arzela theorem, we need to demonstrate that $\{z_i^n (t, x), n\geq 1\}$ is equicontinuous in $C([0, T], L^2(\Omega))$. This is indeed true: because of the boundedness of $\dfrac{\partial z_i^n}{\partial t}$ in $L^2(Q)$, there exists a positive constant $k$ such that

for all $s, t\in [0, T]$. Hence, $z_i^n$ is compact in $C([0, T], L^2(\Omega))$ and there exists a subsequence of $\{z_i^n\}$, denoted also $\{z_i^n\}$, converging uniformly to $z^*_i$ in $L^2(\Omega)$ with respect to $t$. Since $\Delta z^n_i$ is bounded in $L^2(Q)$, there exists a sub-sequence, denoted again $\Delta z^n_i$, converging weakly in $L^2(Q)$. For every distribution $\varphi$,

Thus, $\Delta z^n_i \rightharpoonup \Delta z^*_i$ in $L^2(Q)$. By the same argument, $\dfrac{\partial z_i^n}{\partial t} \rightharpoonup \frac{\partial z^*_i}{\partial t}$ and $z_i^n \rightharpoonup z_i^*$ in $L^2(0, T;H^2(\Omega))$ and $z_i^n \rightharpoonup z_i^*$ in $L^\infty(0, T;H^1(\Omega))$. From $z_1^nz_2^n = (z_1^n-z_1^*)z_2^n+z_1^n(z_2^n-z_2^*)$, we deduce that $z_1^nz_2^n \rightharpoonup z_1^*z_2^*$ in $L^2(Q)$. Therefore, $u_n \rightharpoonup u^*$ in $L^2(Q)$. Since $U_{ad}$ is closed, then $u^* \in U_{ad}$.

Step 3: We conclude that $u^n z_2^n \rightharpoonup u^*z_2^*$ in $L^2(Q)$. Letting $n\rightarrow \infty$ in (5.2), we obtain that $z^*$ is a solution of equation (4.1) corresponding to $u^*$. Therefore,

This shows that $J$ attains its minimum at $(z^*, u^*)$.

6.   Necessary optimality conditions

Now we characterize the optimality that we proved to exist in Section 5. Let $(z^*, u^*)$ be an optimal pair and $u^\epsilon = u^*+\epsilon u$, $\epsilon > 0$, be a control function such that $u \in L^2(Q)$ and $u \in U_{ad}$. We denote by $z^\epsilon = (z_1^\epsilon, z_2^\epsilon, z_3^\epsilon, z_4^\epsilon)$ and $z^* = (z_1^*, z_2^*, z_3^*, z_4^*)$ the corresponding trajectories associated with the controls $u^\epsilon$ and $u^*$, respectively.

In the following result we decompose the right-hand side of our control system into three quantities: $M$, related to the Laplacian part; $R$, linked to the control part; and $F$ for the remaining terms.

Theorem 3. For all $i\in\{1, 2, 3, 4\}$, the mapping $u \longmapsto z_i(u)$ defined from $U_{ad}$ to $W^{1, 2}([0, T], H(\Omega))$ is Gateaux differentiable with respect to $u^*$. For all $u \in U_{ad}$, set $z_i'(u^*)u = Z_i$. Then $Z = (Z_1, Z_2, Z_3, Z_4)$ is the unique solution of the problem

where

Proof. Put $Z^\varepsilon_i = \frac{z^\varepsilon_i-z^*_i}{\varepsilon}$. By subtracting the two systems verified by $z^\varepsilon_i$ and $z^*_i$, we get

Consider the semigroup $(\Gamma(t), t \geq 0)$ generated by $M$. Then the solution of this system is given by

Since the elements of the matrix $F^\varepsilon$ are uniformly bounded with respect to $\varepsilon$, according to Grönwall's inequality one has that $Z^\varepsilon_i$ is bounded in $L^2(Q)$. Hence, $z^\varepsilon_i \rightarrow z^*_i$ in $L^2(Q)$. Letting $\varepsilon \rightarrow 0$, we have

Adopting the same technique, we deduce that $Z^\varepsilon_i \rightarrow Z^*_i$ as $\varepsilon \rightarrow 0$.

Let $p = (p_1, p_2, p_3, p_4)$ be the adjoint variable of $Z$ and denote by $F^*$ the adjoint of the Jacobian matrix $F$. We can write the dual system associated to our problem as

where

Lemma 4. Under the hypothesis of Theorem 1, the system (6.1)of adjoint variables admits a unique solution $p \in W^{1, 2}([0, T], H(\Omega))$ with $p_i \in G(T, \Omega)$, $i = 1, 2, 3, 4$.

Proof. The result follows by the change of variables $s = T-t$ so as to apply the same method performed in the proof of Theorem 3.

We are now in a position to obtain a necessary optimality condition for the optimal control $u^*$.

Theorem 5. If $u^*$ is an optimal control and $z^*\in W^{1, 2}([0, T];H(\Omega))$ is its corresponding solution, then

Proof. Let $u^*$ be an optimal control and let $z^*$ be the corresponding optimal state. Set $u^\varepsilon = u^*+\varepsilon u \in U_{ad}$ and let $z^\varepsilon$ be the corresponding state trajectory. We have

Since $\lim\limits_{\varepsilon\rightarrow 0} \dfrac{z_2^\varepsilon-z_2^*}{\varepsilon} = \lim\limits_{\varepsilon\rightarrow 0} \dfrac{z_2(u^*+\varepsilon h)-z_2^*}{\varepsilon} = Z_2$, $\lim\limits_{\varepsilon\rightarrow 0} z_2^\varepsilon = z_2^*$ and $z_2^\varepsilon, z_2^*\in L^\infty(Q)$, then $J$ is Gateaux differentiable with respect to $u^*$ with

If we take $u = v-u^*$, then we obtain

Since

and $U_{ad}$ is convex, then $J'(u^*)(v-u^*)\geq 0$ for all $v\in U_{ad}$, which is equivalent to

Thus, $bu^* = -{R^*} p$ and, consequently, $u^* = \dfrac{z_2^*(p_2-p_1)}{b}$. Since $u^*\in U_{ad}$, we have that (6.2) holds.

Note that Theorem 5 provides a constructive method, giving an explicit expression (6.2) for the optimal control.

7.   Conclusions and future work

We have extended the time deterministic epidemic SICA model due to Silva and Torres ^[12] to spatiotemporal dynamics, which take into account not only the local reaction of appearance of new infected individuals but also the global diffusion occurrence of the other infected individuals. This allows to incorporate an additional amount of arguments into the system. More precisely, firstly we have modeled the spatiotemporal behavior by incorporating the well-known Laplace operator, which has been employed in the literature, in different contexts, to better understand what happens during any possible displacement of different species and individuals. Here, we justify and interpret its use in the context of HIV/AIDS epidemics. Secondly, we have presented an optimal control problem to minimize the number of infected individuals through a suitable cost functional. Proved results include: existence and uniqueness of a strong global solution to the system, obtained using some adapted tools from semigroup theory; some characteristics of the existing solution; existence of an optimal control, investigated using an effective method based on some properties within the weak topology; and necessary optimality conditions to quantify explicitly the optimal control.

As future work, we plan to develop numerical methods for spatiotemporal optimal control problems, implementing the necessary optimality conditions we have proved here. This is under investigation and will be addressed elsewhere. Another interesting line of research concerns the bifurcation analysis for different parameters.

Acknowledgments

This research was funded by The Portuguese Foundation for Science and Technology (FCT-Fundacão para a Ciência e a Tecnologia), grant number UIDB/04106/2020 (CIDMA). The authors are very grateful to three anonymous Reviewers for several constructive questions and remarks that helped them to improve their work.

Conflict of interest

The authors declare that there are no conflicts of interest.

Appendix

A.   Simulation method and code

The focus of our work is more theoretical, linked to the proposed spatiotemporal SICA epidemic model (3.2). In Section 5, to motivate our study on optimal control, we have incorporated some selected control values in order to present some adequate scenarios showing the dynamic evolution of the system. In our simulations, we have adopted the first order explicit Euler method to discretize the temporal derivatives and the second order explicit Euler method to discretize the Laplacian operator. Follows our $\textsf{Octave}$/$\textsf{Matlab}$ code:

The reader interested in the scientific computing tool $\textsf{GNU Octave}$ or $\textsf{Matlab}$ is referred to ^[7].