1. Introduction

It is worth to start by quoting Bradley [15] "For real progress, the mathematical modeller, as well as the epidemiologist must have mud on his boots"!

Indeed most of the pioneers in mathematical epidemiology have got "mud on their boots"; it is a duty and a pleasure to acknowledge here the ones who, apart from D. Bernoulli (1760) [14], established the roots of this field of research (in chronological order): W. Farr (1840) [41], W.H. Hamer (1906) [45], J. Brownlee (1911) [17], R. Ross (1911) [63], E. Martini (1921) [59], A. J. Lotka (1923) [57], W.O. Kermack and A. G. McKendrick (1927) [51], H. E. Soper (1929) [66], L. J. Reed and W. H. Frost (1930) [42], [1], M. Puma (1939) [62], E. B. Wilson and J. Worcester (1945) [69], M. S. Bartlett (1949) [12], G. MacDonald (1950) [58], N.T.J. Bailey (1950) [11], before many others; the pioneer work by En'ko (1989) [39] suffered from being written in Russian; historical accounts of epidemic theory can be found in [64], [35], [36]. After the late '70's there has been an explosion of interest in mathematical epidemiology, also thanks to the establishment of a number of new journals dedicated to mathematical biology. The above mentioned pioneers explored possible models to match real data, based on genuine epidemiological reasoning; further they did not choose a priori deterministic models as opposed to stochastic models. Unfortunately the most recent literature has suffered of a dramatic splitting in both approaches and methods, which has induced criticism among applied epidemiologists. About the relevance of mathematics in Life Sciences, Wilson and Worcester had since long [69] expressed a fundamental statement that we like to share: "Although mathematics is used to develop the logical inferences from known laws, it may also used to investigate the consequences of various assumptions when the laws are not known, .... one of the functions of mathematical and philosophical reasoning is to keep us alive to what may be only possibilities, when the actualities are not yet known".

The scheme of this presentation is the following: in Section 2 a general structure of mathematical models for epidemic systems is presented in the form of compartmental systems; in Paragraph 2.1 the possible derivation of deterministic models is presented as an approximation, for large populations, of stochastic models; in Paragraph 2.2 nonlinear models are discussed as opposed to the standard epidemic models based on the "law of mass" action assumption; in Paragraph 2.3 the concept of field of forces of infection is discussed for structured populations. In Section 3 the particular case of man-environment-man infection is discussed, and, with respect to these models, in Section 4 optimal control problems are presented in the case of boundary feedback. Finally in Section 5 the most important problem of global eradication via regional control is presented.

---

2. Compartmental models

Model reduction for epidemic systems is obtained via the so-called compartmental models. In a compartmental model the total population (relevant to the epidemic process) is divided into a number (usually small) of discrete categories: susceptibles, infected but not yet infective (latent), infective, recovered and immune, without distinguishing different degrees of intensity of infection; possible structures in the relevant population can be superimposed when required (see e.g. Figure 1).

A key problem in modelling the evolution dynamics of infectious diseases is the mathematical representation of the mechanism of transmission of the contagion. The concepts of "force of infection" and "field of forces of infection" (when dealing with structured populations) will be the guideline of this presentation.

We may like to remark here (see also [19]) that this concept is not very far from the medieval idea that infectious diseases were induced into a human being by a flow of bad air ("mal aria" in Italian). On the other hand in quantum field theory any field of forces is due to an exchange of particles: in this case bacteria, viruses, etc., so that the corpuscular and the continuous concepts of field are conceptually unified.

It is of interest to identify the possible structures of the field of forces of infection which depend upon the specific mechanisms of transmission of the disease among different groups. This problem has been raised since the very first models when age and/or space dependence had to be taken into account.

Suppose at first that the population in each compartment does not exhibit any structure (space location, age, etc.). The infection process ($S$ to $I$) is driven by a force of infection ($f.i.$) due to the pathogen material produced by the infective population and available at time $t$

$(f.i.)(t) = [g(I(.))] (t)$

which acts upon each individual in the susceptible class. Thus a typical rate of the infection process is given by the

$(\text{incidence rate})(t) = (f.i.)(t) S(t).$

From this point of view, the so called "law of mass action" simply corresponds to choosing a linear dependence of $g(I)$ upon $I$

$(f.i.)(t) = k I(t).$

The great advantage, from a mathematical point of view, is that the evolution of the epidemic is described (in the space and time homogeneous cases) by systems of ODE 's which contain at most bilinear terms.

Indeed, for several models of this kind it is possible to prove global stability of nontrivial equilibria. A general result in this direction has been proposed in [13] where it has been shown that many bilinear epidemic systems can be expressed in the general form

${{dz}\over{dt}}=diag(z)(e+Az)+b(z)$

where $z(t) \in \mathbb R_+^n$ is the state vector, $e \in \mathbb R_+^n$ is a constant vector, $A$ is an $n \times n$ constant matrix, and $diag (z)$ is the diagonal matrix with diagonal entries $z_i.$ Further

$b(z)=c+Bz$

with $c\in \mathbb R_+^n$ a constant vector, and $B={\left(b_{ij}\right)}_{i, j=1, \ldots, n}$ a real constant matrix such that

$b_{ij}\geq 0,\;\;\;\; i, j=1, \ldots, n; \;\;\;\; b_{ii}=0,\;\;\;\; i=1, \ldots, n.$

Once a strictly positive equilibrium $z^* \in \mathbb R_+^{n*}$ has been somehow identified, the major "tool" in analyzing these systems is the so called Volterra-Goh Lyapunov function [44],

$V(z):=\sum\limits_{i=1}^n w_i\left(z_i-{z_i^*}-{z_i^*}\, \ln {{z_i}\over{z_i^*}}\right), \;\;\;\; z\in \mathbb R^{n*}$

where $w_i>0, i=1, \ldots, n\, \, ,$ are real constants (the weights).

Here we denote by

$\mathbb R_+^{n*}:=\left\{z \in \mathbb R^n \mid z_i>0,\;\;\;\; i=1, \ldots, n\right\},$

and clearly

$V:=\mathbb R_+^{n*}\rightarrow \mathbb R_+.$

A discussion on $V$ as a relative entropy can be found in [23].

---

2.1. Deterministic approximation of stochastic models

Actually for populations of a limited size, the stochastic version is more appropriate; but it is not difficult to show that for sufficiently large populations, the usual deterministic approximation can be gained via suitable laws of large numbers (see e.g. [40]).

The stochastic process modelling an $S I R$ epidemic, which takes into account the rescaling of the force of infection due to the size of the total population, is a multivariate jump Markov process $(S_t, I_t, R_t)_{t \in \mathbb R_+},$ valued in $\mathbb N^3.$

Considering the usual transitions

$S \rightarrow I \rightarrow R,$

by assuming the law of mass action, the only nontrivial transition rates are usually taken as

$q_{(S, I), (S-1, I+1)} = \kappa \frac{I}{N} S: \;\; \text{infection};$

(1)

$q_{(S, I), (S, I-1)} = \delta I: \;\; \text{removal},$

(2)

$N=S_t+I_t+R_t= S_0 + I_0 = const.$

(3)

We may notice that the above transition rates can be rewritten as follows

$q_{(S, I), (S-1, I+1)} = N \kappa \frac{I}{N} \frac{S}{N};$

(4)

$q_{(S, I), (S, I-1)} = N \delta \frac{I}{N}.$

(5)

So that both transition rates are of the form

$q_{k, k+l}^{(N)}=N \beta_l\left(\frac{k}{N}\right)$

(6)

for

$k=(S, I)$

(7)

and

$k+l=\left\{ \begin{array}{ll} (S, I-1), \\ (S-1, I+1). \end{array}\right.$

(8)

Due to the constancy of the total population we may reduce the analysis to the Markov process $\hat{X}^{(N)}:= (S_t, I_t)$, which satisfies a stochastic evolution equation of the form

$\hat{X}^{(N)}(t)=\hat{X}^{(N)}(0)+\sum\limits_{l\in \mathbb{Z}^2}lY_l\left(N\int _0^t \beta_l\left(\frac {\hat{X}^{(N)}(\tau)}{N}\right)d\tau\right),$

(9)

for $t<\tau_{\infty},$ the possible Markov time of explosion of the epidemic.

Here the $Y_l$ are independent standard Poisson processes, and the sum is carried out only on the $l$'s for which $\beta_l \neq 0.$

By setting

$F(x) =\sum\limits_{l \in \mathbb{Z}^2} l \beta_l(x),\;\;\;\; x \in \mathbb{R}^2$

(10)

for the scaled process

$X^{(N)}=\frac {1}{N} \hat{X}^{(N)},$

(11)

we have

$\begin{align} X^{(N)}(t)=X^{(N)}\, (0) &+ \int _0^t F(X^{(N)}(\tau))d\tau \\ &+ \sum\limits_{l\in \mathbb{Z}^2} \frac {l}{N}\tilde{Y}_l\left(N\int _0^t \beta_l\left(X^{(N)}(\tau)\right)d\tau\right) \end{align}$

(12)

where the

$\tilde{Y}_l(u)=Y_l(u)-u$

(13)

are independent centered standard Poisson processes, so that the last term in the above equation is a zero-mean martingale.

Of interest is the asymptotic behavior of the system for a large value of the scale parameter $N.$

By the strong law of large numbers for Poisson processes (more generally for martingales), we know that

$\lim\limits_{N\to\infty} \sup\limits_{u\leq v}\left|\frac{1}{N}\tilde{Y}_l(Nu)\right|=0, \;\;\;\;\textrm{a.s.},$

(14)

for any $v\geq0$. As a consequence, it is not a surprise the following result, based on Doob's inequality for martingales [40].

Theorem 2.1. Under suitable regularity assumptions on $\beta_l$ and on $F,$ if

$\lim\limits _{N\rightarrow \infty}X^{(N)}(0)=x_0 \in \mathbb{R}^2,$

(15)

then, for every $t\geq 0$,

$\lim\limits _{N\rightarrow \infty} \sup\limits_{\tau\leq t}\left|X^{(N)}(\tau)- x(\tau)\right| =0, \;\;\;\;a.s.,$

(16)

where $x(t)$, $t \in \mathbb{R}_+$ is the unique solution of

$x(t)= x_0 + \int _0^t F(x(s))ds, \;\;\;\;t\geq 0,$

(17)

wherever it exists.

In our case the above deterministic system becomes the usual deterministic SIR model

$\left\{ \begin{array}{ll} \displaystyle{ \frac{d s(t)}{d t}}= - \kappa s(t) s(t) \\ \displaystyle{\frac{d i(t)}{d t}}= \kappa s(t) i(t) - \delta i(t) \end{array} \right.$

(18)

for

$s(t):= \lim\limits_{N \rightarrow \infty} \frac{S_t}{N}, \;\;\;\; i:= \lim\limits_{N \rightarrow \infty} \frac{I_t}{N}.$

A different scaling, may give rise to the diffusion approximation of the epidemic system (see [40], [22], and [67], for a variety of applications to Biology and Medicine).

An interesting "pathology" arises when the relevant populations are very small, so that the deterministic approximation of the epidemic system may fail. Indeed for many epidemic models, above threshold the infective fraction of the relevant deterministic equations, while tending eventually to large values of a possible endemic level, may get very close to zero, but still never becomes extinct. This situation had been analyzed in [47] by suitable perturbation methods on the Fokker-Planck equation associated with the diffusion approximation of a typical SIR epidemic model, which lead to a non trivial extinction probability of the infective population, whenever its deterministic counterpart may get close to zero.

It is worth mentioning that the discussion regarding the original stochastic model and its deterministic counterpart had involved J.L. Doob and others, who proposed (in 1945) [38] an algorithm for generating statistically correct trajectories of the stochastic system. It was presented by D. Gillespie in 1976 [43] as the Doob-Gillespie algorithm, well known in computational chemistry and physics.

---

2.2. Nonlinear models

Referring to the "Law of mass Action", Wilson and Worcester [69] stated the following:

"It would in fact be remarkable, in a situation so complex as that of the passage of an epidemic over a community, if any simple law adequately represented the phenomenon in detail... even to assume that the new case rate should be set equal to any function... might be questioned?".

Indeed Wilson and Worcester [69], and Severo [65] had been among the first epidemic modelers including nonlinear forces of infection of the form

$(f.i.)(t)= \kappa I(t)^p S(t)^q$

in their investigations. Here $I(t)$ denotes the number of persons who are infective, and $S(t)$ denotes the number of persons who are susceptible to the infection.

Independently, during the analysis of data regarding the spread of a cholera epidemic in Southern Italy during 1973, in [28] the authors suggested the need to introduce a nonlinear force of infection in order to explain the specific behavior emerging from the available data.

A more extended analysis for a variety of proposed generalizations of the classical models known as Kermack-McKendrick models, appeared in [29], though nonlinear models became widely accepted in the literature only a decade later, after the paper [55].

Nowadays models with nonlinear forces of infection are analyzed within the study of various kinds of diseases; typical expressions include the so called Holling type functional responses (see e.g. [29], [48])

$(f.i.)(t)=g(I(t))\, \, \, ;$

with

$g(I)=\displaystyle{{k\, \, I^p}\over{\alpha + \beta\, \, \displaystyle{ I^q}}}\, \, \, \, ,\;\;\;\; p, q\, \, >0 \, \, \, .$

(19)

Particular cases are

$g(I)=k\, \, I^p\, \, \, \, ,\;\;\;\; p>0$

(20)

For the case $p=q$ we have the behaviors described in Figure 2.

A rather general analysis regarding existence and stability of nontrivial equilibria for model (19) has been carried out in a series of papers [61], [16], [56], [48] (see also [19], and [68]). The particular case $p=q=2$ in model (19) induces a saddle point behavior as analyzed in [26] and [31] (see Section 3 for the case with spatial structure).

Additional shapes of $g(I),$ as proposed in [29] which may decrease for large values of $I,$ may be interpreted as "awareness" effects in the contact rates. Significant contributions to this concept and related epidemiological issues in recent literature can be found in [37].

Further extensions include a nonlinear dependence upon both $I$ and $S$, as discussed in modelling AIDS epidemics (see e.g. [32], [33], and references therein), where the social structure of the host population is analyzed too.

---

2.3. Structured populations

When dealing with populations which exhibit some structure (identified here by a parameter $z$), either discrete (e.g. social groups) or continuous (e.g. space location, age, etc.), the target of the infection process is a specific "subgroup" $z$ in the susceptible class, so that the force of infection has to be evaluated with reference to that specific subgroup. This induces the introduction of a "field of forces of infection" $(f.i.)(z;t)$ such that the incidence rate at time t at the specific "location" $z$ will be given by

$(incidence \, \, \, \, rate)(z;t)=(f.i.)(z;t)\, \, s(z;t).$

When dealing with populations with space structure the relevant quantities are spatial densities, such as $s(z;t)$ and $i(z;t)$, the spatial densities of susceptibles and of infectives respectively, at a point $z$ of the habitat $\Omega$, and at time $t \geq 0$.

The corresponding total populations are given by

$S(t)=\int_\Omega {s(z;t)\, \, dz},\;\;\;\; I(t)=\int_\Omega {i(z;t)\, \, dz }$

In the law of mass action model, if only local interactions are allowed, the field at point $z\in \Omega$ is given by

$(f.i.)(z;t)=k(z)\, \, i(z;t).$

On the other hand if we wish to take distant interactions too into account, as proposed by D.G. Kendall in [50], the field at a point $z \in \Omega$ is given by

$(f.i.)(z;t)=\int_\Omega {k(z, z')\, \, i(z';t)\, \, dz'}\, \, .$

For this case the emergence of travelling waves has been shown in [50] and [9]. The analysis of the diffusion approximation of Kendall's model can be found in [49].

When dealing with populations with an age structure, we may interpret the parameter $z$ as the age-parameter so that the first model above is a model with intracohort interactions while the second one is a model with intercohort interactions (see e.g. [18], and references therein).

A large literature on the subject can be found in [19].

---

3. Spatially structured man-environment-man epidemics

A widely accepted model for the spatial spread of epidemics in an habitat $\Omega$, via the environmental pollution produced by the infective population, e.g. via the excretion of pathogens in the environment, is the following, as proposed in [20], [21] (see also [19], and references therein).

Typical real cases include typhoid fever, malaria, schistosomiasis, cholera, etc. (see e.g. [34], [6]).

$\left\{ \begin{array}{ll} \displaystyle {{\partial u_1}\over{\partial t}}(x, t)=d_1\Delta u_1(x, t) \displaystyle -a_{11}u_1(x, t)+\int_{\Omega }k(x, x') u_2(x', t)dx' \\ \displaystyle {{\partial u_2}\over{\partial t}}(x, t)= -a_{22}u_2(x, t)+g(u_1(x, t)) \end{array} \right.$

(21)

in $\Omega \subset {\mathbb R}^N\, \ (N\geq 1)$, a nonempty bounded domain with a smooth boundary $\partial \Omega$; for $t \in (0, +\infty )$, where $a_{11}\geq 0$, $a_{22}\geq 0, ~d_1>0$ are constants.

● $u_1(x, t)$ denotes the concentration of the pollutant (pathogen material) at a spatial location $x \in \overline{\Omega },$ and a time $t\geq 0;$

● $u_2(x, t)$ denotes the spatial distribution of the infective population.

● The terms $-a_{11}u_1(x, t)$ and $-a_{22}u_2(x, t)$ model natural decays.

● The total susceptible population is assumed to be sufficiently large with respect to the infective population, so that it can be taken as constant.

Environmental pollution is produced by the infective population, so that in the first equation of System (21), the integral term

$\int_{\Omega }k(x, x') u_2(x', t)dx'$

expresses the fact that the pollution produced at any point $x' \in \Omega$ of the habitat is made available at any other point $x \in \Omega$; when dealing with human pollution, this may be due to either malfunctioning of the sewage system, or improper dispersal of sewage in the habitat. Linearity of the above integral operator is just a simplifying option.

Model (21) includes spatial diffusion of the pollutant, due to uncontrolled additional causes of dispersion (with a constant diffusion coefficient to avoid purely technical complications); we assume that the infective population does not diffuse (the case with diffusion would be here a technical simplification). As such, System (21) can be adopted as a good model for the spatial propagation of an infection in agriculture and forests, too.

The above model is part of another important class of epidemics which exhibit a quasimonotone (cooperative) behavior (see [19]). For this class of problems stability of equilibria can be shown by monotone methods, such as the contracting rectangles technique (see [52], [53]).

The local "incidence rate" at point $x \in \overline{\Omega },$ and time $t\geq 0$, is

$(i.r.)(t)=g(u_1(x, t)),$

depending upon the local concentration of the pollutant.

---

3.1. Seasonality.

If we wish to model a large class of fecal-oral transmitted infectious diseases, such as typhoid fever, infectious hepatitis, cholera, etc., we may include the possible seasonal variability of the environmental conditions, and their impact on the habits of the susceptible population, so that the relevant parameters are assumed periodic in time, all with the same period $T\in (0, +\infty )$.

As a purely technical simplification, we may assume that only the incidence rate is periodic, and in particular that it can be expressed as

$(i.r.)(x, t)=h(t, u_1(x, t))=p(t)g(u_1(x, t)),$

where $h$, the functional dependence of the incidence rate upon the concentration of the pollutant, can be chosen as in the time homogeneous case, with possible behaviors as shown in Figure 2.

The explicit time dependence of the incidence rate is given via the function $p(\cdot),$ which is assumed to be a strictly positive, continuous and $T-$periodic function of time; i.e. for any $t \in \mathbb R,$

$p(t)=p(t+T).$

Remark 1. The results can be easily extended to the case in which also $a_{11}$, $a_{22}$ and $k$ are $T-$periodic functions.

In [21] the above model was studied, and sufficient conditions were given for either the asymptotic extinction of an epidemic outbreak, or the existence and stability of an endemic state; while in [27] the periodic case was additionally studied, and sufficient conditions were given for either the asymptotic extinction of an epidemic outbreak, or the existence and stability of a periodic endemic state with the same period of the parameters.

---

3.2. Saddle point behaviour.

The choice of $g$ has a strong influence on the dynamical behavior of system (21). The case in which $g$ is a monotone increasing function with constant concavity has been analyzed in an extensive way (see [19], [24], [26]); concavity leads to the existence (above a parameter threshold) of exactly one nontrivial endemic state and to its global asymptotic stability. In order to better clarify the situation, consider first the spatially homogeneous case (ODE system) associated with system (21); namely

$\left\{ \begin{array}{ll} \displaystyle {{d z_1}\over{d t}}(t)=\displaystyle -a_{11}z_1(t)+ a_{12} z_2(t) \\ \displaystyle {{d u_2}\over{d t}}(t)= -a_{22}z_2(t)+g(z_1(t)) \end{array} \right.$

(22)

In [26] and [60] the bistable case (in which system (22) may admit two nontrivial steady states, one of which is a saddle point in the phase plane) was obtained by assuming that the force of infection, as a function of the concentration of the pollutant, is sigma shaped. In [60] this shape had been obtained as a consequence of the sexual reproductive behavior of the schistosomes. In [26] (see also [25]) the case of fecal-oral transmitted diseases was considered; an interpretation of the sigma shape of the force of infection was proposed to model the response of the immune system to environmental pollution: the probability of infection is negligible at low concentrations of the pollutant, but increases with larger concentrations; it then becomes concave and saturates to some finite level as the concentration of pollutant increases without limit.

Let us now refer to the following simplified form of System (21), where as kernel we have taken $k(x, x')= a_{12}\delta(x-x'),$

$\left\{ \begin{array}{ll} \displaystyle {{\partial u_1}\over{\partial t}}(x, t)=d_1\Delta u_1(x, t) \displaystyle -a_{11}u_1(x, t)+ a_{12} u_2(x, t) \\ \displaystyle {{\partial u_2}\over{\partial t}}(x, t)= -a_{22}u_2(x, t)+g(u_1(x, t)) \end{array} \right.$

(23)

The concavity of $g$ induces concavity of its evolution operator, which, together with the monotonicity induced by the quasi monotonicity of the reaction terms in (23), again imposes uniqueness of the possible nontrivial endemic state. On the other hand, in the case where $g$ is sigma shaped, monotonicity of the solution operator is preserved, but as we have already observed in the ODE case, uniqueness of nontrivial steady states is no longer guaranteed. Furthermore, the saddle point structure of the phase space cannot be easily transferred from the ODE to the PDE case, as discussed in [26], [31]. In [26], homogeneous Neumann boundary conditions were analyzed; in this case nontrivial spatially homogeneous steady states are still possible. But when we deal with homogeneous Dirichlet boundary conditions or general third-type boundary conditions, nontrivial spatially homogeneous steady states are no longer allowed. In [31] this problem was faced in more detail; the steady-state analysis was carried out and the bifurcation pattern of nontrivial solutions to system (23) was determined when subject to homogeneous Dirichlet boundary conditions. When the diffusivity of the pollutant is small, the existence of a narrow bell-shaped steady state was shown, representing very likely a saddle point for the dynamics of (23). Numerical experiments confirm the bistable situation: "small" outbreaks stay localized under this bell-shaped steady state, while "large" epidemics tend to invade the whole habitat.

---

4. Boundary feedback

An interesting problem concerns the case of boundary feedback of the pollutant, which has been proposed in [24], and further analyzed in [30]; an optimal control problem has been later analyzed in [8].

In this case the reservoir of the pollutant generated by the human population is spatially separated from the habitat by a boundary through which the positive feedback occurs. A model of this kind has been proposed as an extension of the ODE model for fecal-oral transmitted infections in Mediterranean coastal regions presented in [28].

For this kind of epidemics the infectious agent is multiplied by the infective human population and then sent to the sea through the sewage system; because of the peculiar eating habits of the population of these regions, the agent may return via some diffusion-transport mechanism to any point of the habitat, where the infection process is restarted.

The mathematical model is based on the following system of evolution equations:

$\left\{ \begin{array}{ll} \displaystyle{\frac{\partial u_1}{\partial t} }(x;t) =\Delta u_1(x;t)-a_{11} u_1(x;t) \\ \displaystyle{\frac{\partial u_2}{\partial t}}(x;t) =-a_{22}u_2(x;t)+ g(u_1(x;t)) \end{array} \right.$

in $\Omega \times (0, +\infty)$, subject to the following boundary condition

$\displaystyle{\frac{\partial u_1}{\partial \nu}} (x;t) +\alpha u_1(x;t)=\int_{\Omega} k (x, x') u_2(x';t)\, dx'$

on $\partial \Omega \times (0, +\infty)$, and also subject to suitable initial conditions.

Here $\Delta$ is the usual Laplace operator modelling the random dispersal of the infectious agent in the habitat; the human infective population is supposed not to diffuse. As usual $a_{11}$ and $a_{22}$ are positive constants. In the boundary condition the left hand side is the general boundary operator $B:=\displaystyle{{\partial\over{\partial \nu}}+\alpha(\cdot)}$ associated with the Laplace operator; on the right hand side the integral operator

$H\left[u_2(\cdot, t)\right](x):= \int _{\Omega} k(x, x') u_2(x';t)\, dx'$

describes boundary feedback mechanisms, according to which the infectious agent produced by the human infective population at time $t>0$, at any point $x' \in \Omega$, is available, via the transfer kernel $k(x, x')$, at a point $x \in \partial \Omega$.

Clearly the boundary $\partial \Omega$ of the habitat $\Omega$ can be divided into two disjoint parts: the sea shore $\Gamma_1$ through which the feedback mechanism may occur, and $\Gamma_2$ the boundary on the land, at which we may assume complete isolation.

The parameter $\alpha (x)$ denotes the rate at which the infectious agent is wasted away from the habitat into the sea along the sea shore. Thus one may well assume that

$\alpha(x), k(x, \cdot) =0,\;\;\;\; \hbox{ for} \, \, \, \, x\in \Gamma_2\, \, \, .$

A relevant assumption, of great importance in the control problems that we have been facing later, is that the habitat $\Omega$ is "epidemiologically" connected to its boundary by requesting that

$\hbox{ for any} \, \, \, \, \, x'\in\Omega \, \, \, \, \, \hbox{ there exists some}\, \, \, \, \, x\in \Gamma_1 \, \, \, \, \, \hbox{ such that} \, \, k(x, x')>0.$

This means that from any point of the habitat infective individuals contribute to polluting at least some point on the boundary (the sea shore).

In the above model delays had been neglected and the feedback process had been considered to be linear; various extensions have been considered in subsequent literature.

---

5. Regional control: Think Globally, Act Locally

Let us now go back to System (21) in $\Omega \subset {\mathbb R}^N\, \ (N\geq 1)$, a nonempty bounded domain with a smooth boundary $\partial \Omega$; for $t \in (0, +\infty )$, where $a_{11}\geq 0$, $a_{22}\geq 0, ~d_1>0$ are constants.

The public health concern consists of providing methods for the eradication of the disease in the relevant population, as fast as possible. On the other hand, very often the entire domain $\Omega,$ of interest for the epidemic, is either unknown, or difficult to manage for an affordable implementation of suitable environmental sanitation programmes. Think of malaria, schistosomiasis, and alike, in Africa, Asia, etc.

This has led the first author, in a discussion with Jacques Louis Lions in 1989, to suggest that it might be sufficient to implement such programmes only in a given subregion $\omega \subset \Omega,$ conveniently chosen so to lead to an effective (exponentially fast) eradication of the epidemic in the whole habitat $\Omega.$ Though, a satisfactory mathematical treatment of this issue has been obtained only few years later in [2]. This practice may have an enormous importance in real cases with respect to both financial and practical affordability. Further, since we propose to act on the elimination of the pollution only, this practice means an additional nontrivial social benefit on the human population, since it would not be limited in his social and alimentation habits.

In this section a review is presented of some results obtained by the authors, during 2002-2012, concerning stabilization (for both the time homogeneous case and the periodic case). Conditions have been provided for the exponential decay of the epidemic in the whole habitat $\Omega,$ based on the elimination of the pollutant in a subregion $\omega \subset \Omega$. The case of homogeneous third type boundary conditions has been considered, including the homogeneous Neumann boundary conditions (to mean complete isolation of the habitat):

$\frac{\partial u_1}{\partial \nu }(x, t)+\alpha u_1(x, t)=0 \;\;\;\; \mbox{ on } \partial {\Omega}\times (0, +\infty),$

where $\alpha \geq 0$ is a constant and $\partial _{\nu }$ denotes the normal derivative.

For the time homogeneous case the following assumptions have been taken:

(H1) $\, \, g: \mathbb R\rightarrow [0, +\infty )$ is a function satisfying

$a)$ $g(x)=0$, for $x\in (-\infty, 0]$,

$b)$ $g$ is Lipschitz continuous and increasing,

$c)$ $g(x) \leq a_{21}x, \ \mbox{ for any } x\in [0, +\infty )$, where $\, \, a_{21}>0 ;$

(H2) $k\in L^{\infty }(\Omega \times \Omega ), \, \, k(x, x')\geq 0 \,$ a.e. in $\Omega \times \Omega$,

$\int_{\Omega }k(x, x')dx >0 \;\;\;\;\mbox{ a.e. } x'\in \Omega ;$

(H3) $u_1^0, \ u_2^0 \in L^{\infty }(\Omega), u_1^0(x), \ u_2^0(x)\geq0$ a.e. in $\Omega$.

Let $\omega\subset \subset \Omega$ be a nonempty subdomain with a smooth boundary and $\Omega \setminus \overline{\omega }$ a domain. Denote by $\chi_{\omega}$ the characteristic function of $\omega$ (we use the convention

$\chi_{\omega}(x)h(x)=0, \;\;\;\;x\in {\mathbb R}^N \setminus \overline{\omega },$

even if function $h$ is not defined on the whole set ${\mathbb R}^N \setminus \overline{\omega }).$

Our goal is to study the controlled system

$\left\{ \begin{array}{ll} \displaystyle {{\partial u_1}\over {\partial t}}(x, t)= \displaystyle d_1\Delta u_1(x, t)-a_{11}u_1(x, t) +\int_{\Omega }&k(x, x') u_2(x', t)dx' \\ +\chi_{\omega}(x)\, v(x, t), &(x, t)\in \Omega \times (0, +\infty ) \\ \displaystyle {{\partial u_1}\over {\partial \nu }}(x, t)+\alpha u_1(x, t)=0, &(x, t)\in \partial \Omega \times (0, +\infty ) \\ \displaystyle {{\partial u_2}\over {\partial t}}(x, t)=-a_{22}u_2(x, t) +g(u_1(x, t)), \ \ &(x, t)\in \Omega \times (0, +\infty ) \\ u_1(x, 0)=u_1^0(x), u_2(x, 0)=u_2^0(x), &x\in \Omega, \end{array} \right.$

subject to a control $v\in L^{\infty }_{loc}(\overline{\omega }\times [0, +\infty ))$ (which implies that $supp (v(t))\subset {\overline \omega}$ for $t\geq 0$).

We have to mention that existence, uniqueness and nonnegativity of a solution to the above system can be proved as in [10]. The nonnegativity of $u_1$ and $u_2$ is a natural requirement due to the biological significance of $u_1$ and $u_2$.

Definition 5.1. We say that our system is zero-stabilizable if for any $u_1^0$ and $u_2^0$ satisfying (H3) a control $v\in L^{\infty }_{loc}(\overline{\omega }\times [0, +\infty ))$ exists such that the solution $(u_1, u_2)$ satisfies

$u_1(x, t)\geq 0, u_2(x, t)\geq 0, \;\;\;\;\mbox{ a.e. } x\in \Omega, \ \mbox{ for any } t\geq 0$

and

$\displaystyle \lim\limits_{t\to \infty} \Vert u_1(t)\Vert_{L^{\infty }(\Omega)}= \lim\limits_{t\to \infty} \Vert u_2(t)\Vert_{L^{\infty }(\Omega)}=0 .$

Definition 5.2. We say that our system is locally zero-stabilizable if there exists $r_0>0$ such that for any $u_1^0$ and $u_2^0$ satisfying (H3) and $\|u_1^0\|_{L^{\infty }(\Omega )}, \ \|u_2^0\|_{L^{\infty }(\Omega )}\leq r_0,$ there exists $v\in L^{\infty }_{loc}(\overline{\omega }\times[0, +\infty ))$ such that the solution $(u_1, u_2)$ satisfies $u_1(x, t)\geq 0, \ u_2(x, t)\geq 0$ a.e. $x\in \Omega, \ \mbox{ for any } t\geq 0$ and $\lim_{t\rightarrow +\infty}\| u_1(t)\| _{L^{\infty }(\Omega )}= \lim_{t\rightarrow +\infty}\| u_2(t)\| _{L^{\infty }(\Omega )}=0 .$

Remark 2. It is obvious that if a system is zero-stabilizable, then it is also locally zero-stabilizable.

A stabilization result for our system, in the case of time independent $g$, had been obtained in [2]. In case of stabilizability a complicated stabilizing control had been provided. A stronger result (which indicates also a simpler stabilizing control) has been established in [3] using a different approach. Later, in [4] the authors have further extended the main results to the case of a time $T-$periodic function $g$ and provided a very simple stabilizing feedback control.

In [3], by Krein-Rutman Theorem, it has been shown that

$\left\{ \begin{array}{ll} - \displaystyle d_1\Delta \varphi +a_{11}\varphi -{a_{21}\over{a_{22}}}\int_{\Omega }k(x, x') \varphi (x')dx' =\lambda \varphi, \ &x\in \Omega \setminus \overline{\omega }\\ \varphi (x)=0, &x\in \partial \omega \\ \displaystyle {{\partial \varphi }\over {\partial \nu }}(x)+\alpha \varphi (x)=0, &x\in \partial \Omega, \end{array} \right.$

admits a principal (real) eigenvalue $\lambda_{1 }(\omega),$ and a corresponding strictly positive eigenvector $\varphi \in Int (\mathcal K)$ where

${\mathcal K}=\{ \varphi \in L^{\infty }(\Omega); \ \varphi (x)\geq 0 \ \mbox{ a.e. in } \Omega\} .$

The following theorem holds [3]:

Theorem 5.3. If $\lambda _1(\omega )>0$, then for $\gamma \geq 0$ large enough, the feedback control $v:=-\gamma u_1$ stabilizes our system to zero.

Conversely, if $h$ is differentiable at $0$ and $h'(0)=a_{21}$ and if our system is zero-stabilizable, then $\lambda_1(\omega )\geq 0$.

Moreover, the proof of the main result in [3] shows that for a given affordable sanitation effort $\gamma$, the epidemic process can be diminished exponentially if $\lambda _{1, \gamma }^{\omega }>0$, where $\lambda _{1, \gamma }^{\omega }$ is the principal eigenvalue to the following problem:

$\left\{ \begin{array}{ll} - \displaystyle d_1\Delta \varphi +a_{11}\varphi -{a_{21}\over{a_{22}}}\int_{\Omega }k(x, x') \varphi (x')dx' +\gamma \chi _{\omega }\varphi =\lambda \varphi, \ &x\in \Omega \\ \displaystyle {{\partial \varphi }\over {\partial \nu }}(x)+\alpha \varphi (x)=0, &x\in \partial \Omega . \end{array} \right.$

(24)

A natural question related to the practical implementation of the sanitation policy is the following: "For a given sanitation effort $\gamma >0$ in the region $\omega$, is the principal eigenvalue $\lambda _{1, \gamma }^{\omega }$ positive (and consequently can our epidemic system be stabilized to zero by the feedback control $v:=-\gamma u_1$)?"

So, the first problem to be treated is the estimation of $\lambda _{1, \gamma }^{\omega }$. Since this eigenvalue problem is related to a non-self adjoint operator, we cannot use a variational principle (as Rayleigh's for selfadjoint operators); hence in [5] the authors have proposed an alternative method based on the following result:

$\lim\limits_{t\rightarrow +\infty }\int_{\Omega }y^{\omega }(x, t)dx=\zeta -\lambda _{1, \gamma }^{\omega },$

(25)

where $y^{\omega }$ is the unique positive solution to

$\left\{ \begin{array}{ll} \displaystyle {{\partial y}\over {\partial t}}- d_1\Delta y+a_{11}y+\gamma \chi_{\omega }y-{{a_{21}}\over {a_{22}}} \int_{\Omega }k(x, x') y(x', t)dx'&\\ ~~~~~~~~~~~~~~~~~~\displaystyle -\zeta y+(\int_{\Omega }y(x, t)dx)y=0, &x\in \Omega, \ t>0 \\ \displaystyle {{\partial y}\over {\partial \nu }}(x, t)+\alpha y(x, t)=0, &x\in \partial \Omega, \ t>0 \\ y(x, 0)=1, &x\in \Omega, \end{array} \right.$

(26)

and $\zeta >\lambda _{1, \gamma }^{\omega }$ is a constant.

Remark 3. Problem (26) is a logistic model for the population dynamics with diffusion and migration. Since the solutions to the logistic models rapidly stabilize, this means that (25) gives an efficient method to approximate $\lambda _{1, \gamma }^{\omega }$. Namely, for $T>0$ large enough,

$\zeta -\int_{\Omega }y^{\omega }(x, T)dx$

gives a very good approximation of $\lambda _{1, \gamma }^{\omega }$. The above result leads to a concrete numerical estimation of $\lambda _{1, \gamma }^{\omega }$ by analyzing the "large-time" behavior of the system for different values of $\zeta.$

We may also remark that, if in (26)

$y(x, 0)=y_0, \;\;\;\;x\in \Omega,$

(27)

with $y_0$ an arbitrary positive constant, then

$\lim\limits_{t\rightarrow +\infty }\int_{\Omega }y_1^{\omega }(x, t)dx=\zeta -\lambda _{1, \gamma }^{\omega },$

where $y_1^{\omega }$ is the solution to (26)-(27).

Assume now that for a given sanitation effort $\gamma$, the principal eigenvalue to (24) satisfies $\lambda _{1, \gamma }^{\omega }>0$, and consequently $v:=-\gamma u_1$ stabilizes to zero the solution to (21).

Let $\omega _0$ be a nonempty open subset of $\Omega$, with a smooth boundary and such that $\omega _0\subset \subset \Omega$ and $\Omega \setminus \overline{\omega }_0$ is a domain. Consider ${\mathcal O}$ the set of all translations $\omega$ of $\omega _0$, satisfying $\omega \subset \subset \Omega$. Since, after all, our initial goal was to eradicate the epidemics, we are led to the natural problem of "Finding the translation $\omega ^*$ of $\omega$ ($\omega \in {\mathcal O}$) which gives a small value (possibly minimal) of

$R^{\omega }=\int_{\omega }[u_1^{\omega }(x, T)+u_2^{\omega }(x, T)]dx,$

at some given finite time $T>0.$ "

Here $(u_1^{\omega}, u_2^{\omega})$ is the solution of (1.1) corresponding to $v:=-\gamma u_1$, i.e. $(u_1^{\omega}, u_2^{\omega})$ is the solution to

$\left\{ \begin{array}{ll} \displaystyle {{\partial u_1}\over {\partial t}}(x, t)= \displaystyle d_1\Delta u_1(x, t)-a_{11}u_1(x, t) +\int_{\Omega }&k(x, x') u_2(x', t)dx' \\ -\gamma \chi _{\omega }(x) u_1(x, t), &(x, t)\in \Omega \times (0, +\infty ) \\ \displaystyle {{\partial u_1}\over {\partial \nu }}(x, t)+\alpha u_1(x, t)=0, &(x, t)\in \partial \Omega \times (0, +\infty ) \\ \displaystyle {{\partial u_2}\over {\partial t}}(x, t)=-a_{22}u_2(x, t) +g(u_1(x, t)), \ \ &(x, t)\in \Omega \times (0, +\infty ) \\ u_1(x, 0)=u_1^0(x), &x\in \Omega \\ u_2(x, 0)=u_2^0(x), &x\in \Omega . \end{array} \right.$

(28)

For this reason we are going to evaluate the derivative of $R^{\omega }$ with respect to translations of $\omega$. This will allow to derive a conceptual iterative algorithm to improve at each step the position (by translation) of $\omega$ in order to get a smaller value for $R^{\omega }$.

---

5.1. The derivative of $R^{\omega }$ with respect to translations

For any $\omega \in {\mathcal O}$ and $V\in \mathbb R^n$ we define the derivative

$dR^{\omega }(V)=\lim\limits_{\varepsilon \rightarrow 0}{{R^{\varepsilon V+\omega }-R^{\omega }}\over {\varepsilon }}.$

For basic results and methods in the optimal shape design theory we refer to [46].

Theorem 5.4. For any $\omega \in {\mathcal O}$ and $V\in \mathbb R^n$ we have that

$dR^{\omega }(V)=\gamma \int_0^T\int_{\partial \omega }u_1^{\omega }(x, t)p_1^{\omega }(x, t)\nu (x)\cdot Vd\sigma \ dt,$

where $(p_1^{\omega }, p_2^{\omega })$ is the solution to the adjoint problem

$\left\{ \begin{array}{ll} \displaystyle {{\partial p_1}\over {\partial t}}+d_1\Delta p_1-a_{11}p_1-\gamma \chi_{\omega }p_1+g'(u_1^{\omega })p_2=0,&x\in \Omega, \ t>0 \\ \displaystyle {{\partial p_2}\over {\partial t}} + \int_{\Omega }k(x', x)p_1(x', t)dx'-a_{22}p_2 =0,&x\in \Omega, \ t>0 \\ \displaystyle {{\partial p_1}\over {\partial \nu }}(x, t)+\alpha p(x, t)=0, &x\in \partial \Omega, \ t>0 \\ p_1(x, T)=p_2(x, T)=1, &x\in \Omega . \end{array} \right.$

(29)

Here $\nu (x)$ is the normal inward versor at $x\in \partial \omega$ (inward with respect to $\omega$).

For the construction of the adjoint problems in optimal control theory we refer to [54].

Based on Theorem 5.4, in [5] the authors have proposed a conceptual iterative algorithm to improve the position (by translation) of $\omega \in {\mathcal O}$ (in order to obtain a smaller value for $R^{\omega }$.

---

5.2. The periodic case

As a purely technical simplification, we have assumed that only the incidence rate is periodic, and in particular that it can be expressed as

$(i.r.)(x, t)=h(t, u_1(x, t))=p(t)g(u_1(x, t)),$

were $g$, the functional dependence of the incidence rate upon the concentration of the pollutant, can be chosen as in the time homogeneous case.

In this case our goal is to study the controlled system

$\left\{ \begin{array}{ll} \displaystyle{\frac{\partial u_1 }{\partial t }}(x, t)= \displaystyle d_1\Delta u_1(x, t)-a_{11}u_1(x, t) +\int_{\Omega } &k(x, x') u_2(x', t)dx'\\ + \chi_{\omega}(x) v(x, t), &(x, t)\in \Omega \times (0, +\infty ) \\ \displaystyle{\frac{\partial u_1 }{\partial \nu }}(x, t)+\alpha u_1(x, t)=0, &(x, t)\in \partial \Omega \times (0, +\infty ) \\ \displaystyle{\frac{\partial u_2 }{\partial t }}(x, t)=-a_{22}u_2(x, t) +h(t, u_1(x, t)), \ \ &(x, t)\in \Omega \times (0, +\infty ) \\ u_1(x, 0)=u_1^0(x), u_2(x, 0)=u_2^0(x), &x\in \Omega, \end{array} \right.$

(30)

with a control $v\in L^{\infty }_{loc}(\overline{\omega }\times [0, +\infty ))$ (which implies that $supp (v(\cdot, t))\subset {\overline \omega}$ for $t\geq 0$).

The explicit time dependence of the incidence rate is given via the function $p(\cdot),$ which is assumed to be a strictly positive, continuous and $T-$periodic function of time; i.e. for any $t \in \mathbb R,$

$p(t)=p(t+T).$

Remark 4. The results can be easily extended to the case in which also $a_{11}$, $a_{22}$ and $k$ are $T-$periodic functions.

Consider the following (linear) eigenvalue problem

$\left\{ \begin{array}{ll} \displaystyle{\frac{\partial \varphi }{\partial t}}- \displaystyle d_1\Delta \varphi +a_{11}\varphi-\int_{\Omega }k(x, x') \psi(x', t)dx' =\lambda \varphi, &x\in \Omega \setminus \overline{\omega }, \ t>0\\ \displaystyle{\frac{\partial \varphi }{\partial \nu}} (x, t)+\alpha \varphi (x, t)=0, &x\in \partial \Omega, \ t>0\\ \varphi (x, t)=0, &x\in \partial \omega, \ t>0\\ \displaystyle{\frac{\partial \psi }{\partial t}} (x, t)+a_{22}\psi (x, t)-a_{21}p(t)\varphi (x, t)=0, &x\in \Omega \setminus \overline{\omega }, \ t>0\\ \varphi (x, t)=\varphi (x, t+T), \psi (x, t)=\psi (x, t+T), &x\in \Omega \setminus \overline{\omega }, \ t\geq 0. \end{array} \right.$

(31)

By similar procedures, as in the time homogeneous case, Problem (31) admits a principal (real) eigenvalue $\lambda _{1 }^T(\omega),$ and a corresponding strictly positive eigenvector $\varphi^T \in Int (\mathcal K^T)$ where

${\mathcal K^T}=\{ \varphi \in L^{\infty }(\Omega \times (0, T)); \ \varphi (x, t)\geq 0 \ \mbox{ a.e. in } \Omega \times (0, T)\} .$

Theorem 5.5. If $\lambda _1^T(\omega )>0$, then for $\gamma \geq 0$ large enough, the feedback control $v:=-\gamma u_1$ stabilizes (30) to zero.

Conversely, if $g$ is differentiable at $0$ and $g'(0)=a_{21},$ and if (30) is zero-stabilizable, then $\lambda _1^T(\omega )\geq 0$.

Theorem 5.6. Assume that $g$ is differentiable at $0$. Denote by $\tilde{\lambda }_1^T(\omega )$ the principal eigenvalue of the problem

$\left\{ \begin{array}{ll} \displaystyle{\frac{\partial \varphi }{\partial t}}- \displaystyle d_1\Delta \varphi +a_{11}\varphi -\int_{\Omega }k(x, x') \psi(x', t)dx' =\lambda \varphi, &x\in \Omega \setminus \overline{\omega }, \ t>0\\ \displaystyle{\frac{\partial \varphi }{\partial \nu}} (x, t)+\alpha \varphi (x, t)=0, &x\in \partial \Omega, \ t>0\\ \varphi (x, t)=0, &x\in \partial \omega, \ t>0\\ \displaystyle{\frac{\partial \psi }{\partial t}} (x, t)+a_{22}\psi (x, t)-g'(0)p(t)\varphi (x, t)=0, &x\in \Omega \setminus \overline{\omega }, \ t>0\\ \varphi (x, t)=\varphi (x, t+T), \psi (x, t)=\psi (x, t+T), &x\in \Omega \setminus \overline{\omega }, \ t\geq 0 \end{array} \right.$

(32)

If $\tilde{\lambda }_1^T(\omega )>0$, then the system is locally zero stabilizable, and for $\gamma \geq 0$ sufficiently large, $v:=-\gamma u_1$ is a stabilizing feedback control.

Conversely, if the system is locally zero stabilizable, then $\tilde{\lambda }_1^T(\omega )\geq 0$.

Remark 5. Since $g'(0)\leq a_{21}$, it follows that $\lambda _1^T(\omega ) \leq \tilde{\lambda }_1^T(\omega ).$ We conclude now that

$1^0$ If $\lambda _1^T(\omega )>0$, the system is zero-stabilizable;

$2^0$ If $\tilde{\lambda }_1^T(\omega )>0$ and $\lambda _1^T(\omega )\leq 0$, the system is locally zero-stabilizable;

$3^0$ If $\tilde{\lambda }_1^T(\omega )<0$, the system is not locally zero-stabilizable and consequently it is not zero stabilizable.

Remark 6. Future directions. Another interesting problem is that when $\omega$ consists of a finite number of mutually disjoint subdomains. The goal is to find the best position for each subdomain. A similar approach can be used.

In a recently submitted paper [7], the problem of the best choice of the subregion $\omega$ has been faced for a general harvesting problem in population dynamics as a shape optimization problem; our future aim is to apply those results to our problem of eradication of spatially structured epidemics.

---

Acknowledgments

The work of V. Capasso was supported by the MIUR-PRIN grant $200777BWEP-003$: "From the stochastic modelling to the statistics of space-time structured population dynamics". The work of Sebastian Aniţa was supported by the CNCS-UEFISCDI (Romanian National Authority for Scientific Research) grant 68/2.09.2013, PN-II-ID-PCE-2012-4-0270: "Optimal Control and Stabilization of Nonlinear Parabolic Systems with State Constraints. Applications in Life Sciences and Economics".

It is a pleasure to acknowledge the contribution by Klaus Dietz regarding the bibliography on the historical remarks reported in the Introduction.

Thanks are due to the Anonymous Referees for their precious advise and suggestions.

---