Environmental health risk analysis from detergent contamination in well water

1.   Introduction

Vector-borne diseases exhibit a common feature: a two-way transmission of the pathogen between human hosts and various species of arthropods, most often mosquitos, that serve as carriers (or vectors) of the pathogen, but direct host-host transmission is impossible for most of them. Diseases caused by pathogens transmitted by mosquito bites (malaria, dengue fever, yellow fever, zika, chikungunya) cause periodic outbreaks in the countries located in tropical climate zones, and pose a growing burden on their healthcare systems, as seen, for instance, in a 30-fold increase in dengue fever's incidence since the 1970s ^[1]. Control campaigns against them are based on two broad types of measures: suppression of transmission and elimination of the vector. They rely on the use of long-lasting insecticidal nets (LLINs), indoor residual spraying with insecticide, use of larvicide, as well as use of personal protection means with repellent activity that include spray-on repellents or repellent-treated textiles and household items.

The success of control measures aimed at the vector's demography, like bed nets and insecticide spraying, in reducing the disease burden has been limited for various reasons. On the one hand, there is evidence for changing patterns in the transmission of malaria to outdoor environments ^[2,3,4,5], while the main vectors of dengue, yellow fever and zika viruses, female mosquitoes of the Aedes genus are active during daytime hours. Indiscriminate use of insecticide and larvicide, on the other hand, poses danger to human health and the ecosystem. Furthermore, insecticide resistance has contributed to the failure of the Global Malaria Eradication Programme ^[6] and has been noted among Aedes mosquitos ^[7,8].

Use of repellent-treated clothes and household items as preventive measure could serve to suppress transmission of pathogens from vector to host and vice versa. Changing mosquito host-seeking behavior by repellent application could reduce the biting rate of the female mosquitoes, and this strategy has gained attention among entomologists. This includes studies of the effect of spray-on repellents and wearable repellent devices on mosquito host-seeking behavior ^[9], and tests of fibers for controlled release of volatile mosquito repellents to prevent bites in the context of malaria control ^[10]. Development of mosquito repellents based on organic substances from plant oils ^[11], and of methods for estimation of efficacy of repellents used in textiles and paints ^[12] reveal new perspectives for such measures.

Mathematical modeling can help analyze and predict the performance of such interventions for meeting policy objectives. We study a simple Ross-Macdonald type model ^[6] for a vector-borne disease with control based on distribution of textiles and household items with repellent properties. The constraints on an intervention campaign of this kind encompass, first, the economic cost of production, distribution, etc. of the repellent-treated textiles, which influences its extent, that is, the maximum fraction $\bar u$ of the target population employing products with mosquito repellent property. Second, due to technological and physical limitations, such as evaporation, washing, UV radiation exposure, etc., the repellent property is assumed to be lost after a period of duration $T$ days..

In this study we propose a method to analyze whether a control campaign under those constraints is able to reduce and maintain the size of the infected host compartment below a certain level, denoted henceforth as the infection cap.

The infection cap $\bar I$ is a constraint on a state variable of the model (infected hosts) and its value depends on factors, such as the availability of treatment, healthcare system capacity or societal or political tolerances, etc. For the purpose of our analysis, we state the following Questions.

Q.1 Given $\bar u, \bar I$ and initial data on infected hosts and vectors, establish whether the size of the infected host compartment can remain below the infection cap for all future times and find the optimal strategy to maintain it.

Q.2 If the initial data are such that the objective from Q.1 can not be met, we address the following: given $\bar u, \bar I, T$, find the optimal strategy to bring the size of the infected host compartment below the infection cap $\bar I$ in a minimal time not exceeding $T$, and maintain it below the cap $\bar I$ until the campaign's end.

Both questions concern characterization of initial data inside the model's state space unlike standard control analysis in the context of vector-borne diseases, which addresses optimal resource allocation ^{[13,14,15,16]} These questions can be addressed using the toolbox of optimal control theory.

Question Q.1 is a viability problem and we study the existence of viable trajectories of the dynamical system, namely those that satisfy the constraint stated in Q.1. The maximum set of initial data for Q.1 where this constraint is met, is the viability kernel. There has been some analysis of viable controls in epidemiology: in reference ^[17] viability kernels for a model with control on the mosquito population by fumigation have been computed and validated with data from a dengue outbreak in Cali, Colombia. In a previous work ^[18] we have characterized the viability kernel of a model for a vector-borne disease with susceptible, infected and removed compartments for the host population. Even though this problem has an infinite time horizon, its solution shall be used for the solution of Question Q.2.

For Question Q.2 we study a reachability problem for a target set (the viability kernel), and compute the minimal entry time function to it. Such analysis can serve as a guide to policy makers to determine determine whether and how fast the repellent-based control campaign constrained by its coverage $\bar u$ and duration $T$ is able to reduce the size of the infected host compartment below the infection cap $\bar I$. If the minimal time lies within the given time horizon $(0, T)$, the campaign is considered effective. Else, the campaign should be considered ineffective because it cannot meet this objective. The reachability problem framed by Question Q.2 is important as it can provide an indicator of the control campaign's performance under the constraints to decision makers.

We use a variational approach to compute numerically approximations of the viability kernel of the controlled system and the minimal entry time function to the viability kernel, using sub-zero level set of a value function solving a Hamilton-Jacobi equation associated to the model system. The viability kernel can also be described analytically by its boundary in the phase plane of the state variables (infected hosts and infected vectors) ^[17]. However, the numerical approximation of the value function from the variational approach {allows us} to reconstruct the control function which solves the viability and reachability problems. We illustrate the application of this method from optimal control theory using parameterized models for malaria transmission in Botswana and Zimbabwe ^[19] and discuss the efficacy of the repellent-based control campaign in meeting the stated objectives.

2.   Mathematical model

We use the compartmental model for a vector-borne disease as presented in ^[20], which follows the dynamics of susceptible and infected hosts $S_h(t), I_h(t)$ and susceptible and infected female mosquitoes $S_v(t), I_v(t)$ that serve as vectors for the pathogen, at time $t\ge0$. Susceptible hosts become infected after receiving a bite from an infectious mosquito. A susceptible mosquito becomes infected after biting an infected host, and after an incubation period of length $\tau$ becomes infectious and is able to transmit the pathogen. We model the incubation period without a specific compartment as done by ^[19,21]. If the within-vector incubation period $\tau$ and the expected life-time of a mosquito $1/\mu$ are of the same order of magnitude, the probability of death of an infected mosquito during the pathogen incubation is not negligible. Hence, merely a fraction $\exp(-\mu\tau)$ of the infected mosquitoes $I_v$ is capable of transmitting the pathogen to the host ^[21]. In this model no protective immunity is assumed. Hosts from the infected compartment after a short infectious period re-enter immediately the susceptible compartment ^[22,23], the length of the period is $1/\gamma$, with $\gamma$ the recovery rate.

Let $u(t)$ denote the proportion of the target population at time $t$ using the measures provided by the control campaign. With $k$ being the repellent efficacy, the modified mosquito biting rate due to use of control measures is $a_m(1- k u(t))$. Such modified infection rates are used in reference ^[15] in the context of an age-structured model for malaria and controls by long-lasting insecticidal nets. We assume that the cost of control effort $C(u(t))$ (production and distribution of the repellent-treated textiles or spray-on repellents) is borne by decision makers and is a linear function of $u(t)$ (compare reference ^[38]), and that $C(u(t))\le C_{max}$ for all $t\ge0$. This means that $u(t)\le \bar u$ for all $t\ge0$ for some maximum coverage $\bar u\in(0, 1]$ of the target population. The set of control functions we consider is

Observe that for $u = \bar u$ we have the maximum available reduction of the biting rate, while for $u = 0$ there is no control used, and hence no reduction of the biting rate.

The host population is assumed to be constant over time, $S_h(t)+I_h(t) = N$. We assume that the repellent-based control does not lead to changes in the total populations of female mosquitoes over time, $M = S_v(t)+I_v(t)$. Hence, we can reduce the number of state variables by letting $S_h = N-I_h, S_v = M-I_v$ and non-dimensionalizing $x_1 = {I_h}/{N}, x_2 = {I_v}/{M}$. In our notation, $x(t)$ denotes the vector of state variables describing the different compartments of the epidemiological model at time $t$, and $f(x, u)$ denotes the model dynamics subject to a control function $u\in\mathcal{U}$.

Letting $\nu = {M}/{N}$ be the average number of female mosquitoes per host and setting for shortness's sake $\alpha = a_m p_h\exp(-\mu\tau), \beta = a_m p_v$, we work with the 2-dimensional non-autonomous system:

subject to initial conditions $x^0 = (x_1(0), x_2(0))\in\Omega = [0, 1]^2$. The domain of definition $\Omega$ corresponds to the states $x$ with relevant epidemiological meaning for (2.1). On the segments on the boundary $\partial\Omega$ defined by $x_1 = 0, x_2 = 0, x_2 = 1, x_1 = 1$, one can verify that the flow defined by $f$ points towards $\mathtt{int}\; \Omega$. Thus, $\Omega$ is forward invariant under $f$.

The associated trajectory of (2.1) for a given control function $u\in\mathcal{U}$ with initial condition $x^0$ will be denoted by $\bf{x}^u(t; x^0) = \{(x_1^{u}(t), x_2^{u}(t)), t > 0\}$. The set of feasible trajectories on the time interval $[0, +\infty)$ starting at $x^0\in\Omega$ at time $t = 0$ will be denoted by

The system (2.1) is cooperative (or quasimonotone ^[20]) because $\partial_{ x_i} f_j\ge0, i\neq j$. Hence, tools from the theory of ordinary differential equations such as comparison principles ^[24] are at our disposal to establish properties of the feasible trajectories.

We conclude this section with

Lemma 1. The function $f$ in system (2.1) is Lipschitz continuous on $\Omega$ with Lipschitz constant

Proof. The estimate follows from direct computation of the Euclidean norm for $|f(x, u)-f(x', u')|, \; x, x'\in\Omega, u, u'\in[0, \bar u]$.

Using the Lipschitz continuity of $f$ on $\Omega$ (Lemma 1) as well as the facts that a) the set $\{f(x, u), 0\le u\le\bar u\}$ is convex for all $x\in\Omega$ and b) we can choose $C > 0$ so $f(x, u)\le C(1+|x|)$ on $\Omega$ ($f$ has linear growth), we establish that solutions to system (2.1) exist for $t\ge0$ [25,Chapter III.5]. In addition, since the domain and graph of $f$ are closed, $f$ is a Marchaud map on $\Omega$ [26,Corollary 2.2.5].

3.   Viability analysis

In Question Q.1 we are interested in finding the initial conditions $x^0$ of (2.1) such that for a given maximum extent $\bar u$ and infection cap $\bar I > 0$, there exists $u(t)\in\mathcal{U}$ such that $x_1^u(t)\le \bar I, \forall t\ge0$. Let ${\mathcal A}(\bar I) = \{(x_1, x_2)\in\Omega \; |\; x_1\le \bar I\}$, which is a closed subset of $\Omega$ with a compact boundary. A feasible trajectory $\bf{x}\in{\mathcal S}(x^0, \mathcal{U})$ is called viable for $\mathcal{A}(\bar I)$ if the constraint $x_1^u(t)\le \bar I, \forall t\ge0$ is met, in other words, $\bf{x}(t)\in \mathcal{A}(\bar I)$ for $t \in [0, +\infty]$. A set $\mathcal{D}\subset\Omega$ is a viability domain for $f$ if from any initial state $x^0\in\mathcal{D}$, at least one trajectory $\bf{x}\in{\mathcal S}(x^0, \mathcal{U})$ is viable for $\mathcal{D}$ ^[26,27].

To answer Question Q.1, we characterize the viability kernel associated to $\mathcal{A}\bar I, \bar u$, that is the following set of initial data for (2.1):

Evidently, ${\mathbb V}(\bar I, \bar u)$ is the set of all initial data $x^0\in\mathcal A(\bar I)$ from which viable trajectories for $\mathcal{A}(\bar I)$ for $t\ge0$ exist. It is the largest closed viability domain of $f$ inside $\mathcal{A}(\bar I)$ ^[27]. The set $\mathcal{A}(\hat I)$ being closed, and $f$ being a Marchaud map, the viability kernel ${\mathbb V}(\bar I, \bar u)$ for system (2.1) is well-defined [26,Theorem 4.1.2]. The viability kernel will depend on the constraint on the control $\bar u$, and on the state constraint (the infection cap $\bar I$).

We introduce the effective basic reproduction number for control intervention $u(t)\equiv\bar u$

and the critical threshold of infection

whenever $\mathcal{R}_{\bar u} > 1, I_{crit} > 0$. Following reference ^[18], $I_{crit}$ equals the share of infected hosts at the endemic equilibrium $\mathcal{E}_{\ast}$ of system (3.4) at the maximum coverage with repellent-treated products.

Here we briefly recall stability properties of the model equilibria under constant controls $u(t) = \bar u$. In that case the system (2.1) becomes autonomous,

System (3.4) has a trivial, disease-free equilibrium at the origin ${\mathcal O}$ and an endemic equilibrium:

Note that $\mathcal{E}_{\ast}\in\mathtt{int}\; \Omega$^* if and only if $\mathcal{R}_{\bar u} > 1$. At $\mathcal{R}_{\bar u} = 1$, a transcritical bifurcation occurs and $\mathcal{E}_{\ast}$ and ${\mathcal O}$ coincide ^[23].

^*By $\mathtt{int}\; X$ we denote the interior of a set $X$.

Lemma 2. The system (3.4) has the following asymptotic behavior:

(i) $\mathcal{R}_{\bar u}\le1$ implies that the disease-free equilibrium ${\mathcal O}$ is globally asymptotically stable for (3.4) in $\Omega$.

(ii) $\mathcal{R}_{\bar u} > 1$ implies that the endemic equilibrium $\mathcal{E}_{\ast}$ (3.5) is globally asymptotically stable for (3.4) in $\Omega\setminus{\mathcal O}$.

The proofs of Lemma 2 and further statements are in the Supplementary Material.

Using $\bar I$ and $I_{crit}$, we characterize whether the viability kernel for system (2.1) has a positive Lebesgue measure in $\mathbb{R}^2$:

Proposition 1. Consider the model with constant vector population (2.1). Let

(i) $\bar I \ge \max\{0, I_{crit}\}$. Then $\mathtt{int}\; {\mathbb V}(\bar I, \bar u)$ has positive Lebesgue measure in $\mathbb{R}^2$.

(ii) $I_{crit} > 0$ and $\bar I < I_{crit}$. Then ${\mathbb V}(\bar I, \bar u) = {\mathcal O}$.

Based on Proposition 1, the description of ${\mathbb V}(\bar I, \bar u)$ can be refined by defining its boundaries. This has been the approach in reference ^[17], and we recall the details for model (2.1). We start by a technical result:

Lemma 3. Let $\bar I > max\{0, I_{crit}\}$. Set $\bar v = \frac{\gamma \bar I}{\alpha\nu(1- k \bar u)(1-\bar I)}$ and assume $\bar v < 1$.The initial value problem

has a unique non-negative solution which is monotone decreasing on the interval $(\bar v, \min\{1, v^\ast\})$, where $v^\ast$ is such that $z(v^\ast) = 0$.

Due to the inverse function theorem, Lemma 3 implies the existence of an inverse map $z^{-1}$ to $z$ defined on $[0, \bar I]$ such that $z^{-1}(z(s)) = s$. We denote by $\mathcal{Z}$ the solution curve $(z(s), s)$ in $\Omega$ of the problem in Lemma 3 to describe explicitly the boundary of the viability kernel for the case $\bar I > \max\{0, I_{crit}\}$.

Proposition 2. For given $\bar I, \bar u$ the viability kernel ${\mathbb V}(\bar I, \bar u)$ takes one of the following forms

(i) Let $\bar I > \frac{(1- k \bar u)\alpha\nu}{(1- k \bar u)\alpha\nu+\gamma}$, then ${\mathbb V}(\bar I, \bar u) = \mathcal{A}(\bar I)$.

(ii) Let $\bar I\in(\max\{0, I_{crit}\}, \frac{(1- k \bar u)\alpha\nu}{(1- k \bar u)\alpha\nu+\gamma}]$. Then

where $z(x)$ is the solution to the initial value problem (3.6) in Lemma 3, and $v^\ast$ such that $z(v^\ast) = 0$.

Unfortunately, this "direct" approach does not work for the case $\bar I = I_{crit}$. While Proposition 1 shows that the viability kernel has positive measure, the approach from Lemma 3 cannot work, because the right-hand side of (3.6) is not defined at the initial condition $z(\bar v) = \bar I$, which is precisely at the endemic equilibrium $\mathcal{E}_\ast$.

To avoid this problem, we state an alternative characterization of the viability kernel using a variational formulation. Following reference ^[28], we consider an infinite horizon minimization problem with exact penalization of the state constraint. The constraint is expressed via the cost function, defined by the signed distance function $\Gamma$

Note that $\Gamma(x)\le0$ if and only if $x\in\mathcal{A}(\bar I)$. Following the steps in reference ^[28], we consider for $\ell > 0$ the value function

Note that $v$ takes finite values on $\Omega$. From (3.9), it follows $v(x)\le 0$ if and only if $\bf{x}^u(t; x)\in{\mathbb V}(\bar I, \bar u)$ for all $t\in(0, +\infty)$. Therefore, the viability kernel may be characterized using sub-zero-level sets of the value function $v$:

Let $\ell > L$ with $L$ being the Lipschitz constant of $f$ (Lemma 1). Then the value function $v$ satisfies the following dynamic programming principle

and is Lipschitz continuous ^[28].

Let $\nabla v$ denote the gradient of $v$, $\nabla v = (\partial_{x_1}v, \partial_{x_2}v)$. The results in reference [28,Section 4] imply that $v$ is the unique continuous viscosity solution [33,Definition 2] of the Hamilton-Jacobi-Bellman equation

with $\ell$ as above and Hamiltonian

We refer the reader to reference [25,Chapter 1.3] for the definitions of viscosity solutions to partial differential equations.

Hence, using well-known numerical methods for Hamilton-Jacobi equations based on finite difference discretization of the spatial derivative ^[29,30], one can compute the solution of (3.11) and find an approximation of the viability kernel ${\mathbb V}(\bar I, \bar u)$ for given values of target infection cap $\bar I$ and maximum coverage $\bar u$.

In our case the quantity $\max_{u\in\mathcal{U}}{\mathcal H}(x, u, \nabla v)$ has an explicit form :

Using the computed value function $v$, we can answer Q.1, compute the optimal control function $u(t)$ and reconstruct the optimal trajectory for a given initial condition $x^0 \in {\mathbb V} (\bar I, \bar u)$ using the reconstruction algorithm from reference ^[31,32]. This algorithm is applied in reference ^[18] to approximate the viability kernel for another model of a vector-borne disease.

4.   Minimal entry time problem

Once we have computed the viability kernel, we turn our attention to the reachability problem from Question Q.2. Let the initial data $x^0 = (x_1(0), x_2(0))$ on infected hosts and vectors, the maximum extent $\bar u$, the duration $T$, and the infection cap $\bar I$ be given. We ask if there exists an optimal strategy $u(t)$ to bring $x_1$ population below the infection cap $\bar I$ in a minimal time $\tau < T$ and maintain $x_1(t) < \bar I$ for $t\in(\tau, T]$.

To this purpose we define the minimal entry time function $\mathcal{T}_X(x)$ for initial data $x$ to the target set $X$ as follows:

The objective in Question Q.2 is to compute the minimal entry time function for the target set being the viability kernel ${\mathbb V}(\bar I, \bar u)$. This makes sense unless the viability kernel is trivial (in other words, we are interested in computing $\mathcal{T}$ for ${\mathbb V}(\bar I, \bar u)\neq{\mathcal O}$). By definition, ${\mathbb V}(\bar I, \bar u)\subset\Omega$ is a compact set, hence $\partial{\mathbb V}(\bar I, \bar u)$ is compact. We shall characterize the minimal entry time function using the value function $v$ from (3.9) based on the approach presented in references ^[28,33].

Using the value function $v$ which solves (3.11) and defines the viability kernel ${\mathbb V}(\bar I, \bar u)$ via its sub-zero level set (3.10) from the previous section we formulate a variational problem for the minimal entry time function. The backwards reachable set for the target set (which is viability kernel ${\mathbb V}(\bar I, \bar u)$) at time $T$, denoted by ${\mathbb W}(T)$, is defined as the set of those initial data $x$ for system (2.1) such that there exists a feasible trajectory $\bf{x}\in\mathcal{S}(x, \mathcal{U})$ with $\bf{x}^u(t; x)\in {\mathbb V}(\bar I, \bar u)$ for $t\le T$.

In other words,

Note that $\{{\mathbb W}(t), t\ge0\}$ is a family of increasing closed sets [33,Remark 2]. We remark that even though the biologically relevant domain for system (2.1) is $\Omega$, the proof of Lemma 2 implies $\Omega$ attracts the trajectories of (2.1) starting at any $x^0\in\mathbb{R}_+^2$, so in principle there exist $T > 0$ with $\Omega\subset{\mathbb W}(T)$.

The set ${\mathbb W}(T)$ can be characterized using a level-set approach based on the computed value function $v$ by avoiding any controllability assumption for (2.1). We formulate a minimization problem (4.3) for a value function $w$, which is the unique continuous viscosity solution of a Hamilton-Jacobi-Bellman equation [33,Definition 2 and Theorem 2].

To impose a restriction ${\mathbb W}(T)\subseteq\Omega$ into the variational formulation, we introduce a penalization term $G_\Omega(x)$, which is a Lipschitz-continuous function which satisfies the condition for sub-zero-level set, $G_\Omega(x)\le 0$ if and only if $x\in\Omega$.

To compute the minimal entry time function, we consider the solution $w$ of the minimization problem

Note that $w$ is defined over $\mathbb{R}^2$, but the penalization term in (4.3) makes the value function $w$ positive if the trajectory leaves the domain $\Omega$. The backwards reachable set ${\mathbb W}(T)$ for (2.1) is thus defined as the sub-zero-level set of the function $w$, ${\mathbb W}(T) = \{x\in\mathbb{R}^2\; |\; w(T, x)\le0\},$ and (4.3) implies ${\mathbb W}(T)\subseteq\Omega$.

Following reference [33,Lemma 1], the value function $w$ satisfies the following dynamic programming principle

and $w$ is the unique continuous viscosity solution to the Hamilton-Jacobi-Bellman equation

with $\nabla w = (\partial_{x_1}w, \partial_{x_2}w)$, and Hamiltonian ${\mathcal H}$ provided by (3.12).

Using $w$, we define the minimal entry time function (4.1) $\mathcal{T}_{{\mathbb V}(\bar I, \bar u)}(x)$ for ${\mathbb V}(\bar I, \bar u)$ for a trajectory with initial condition $x$ as

As there is no controllability assumption for (2.1), $\mathcal{T}_{{\mathbb V}(\bar I, \bar u)}$ need not be continuous over $\Omega$. In fact, the proof of part (ii) of Proposition 2 reveals that the minimal entry time function exhibits a discontinuity across the boundary of the viability kernel ${\mathbb V}(\bar I, \bar u)$ traced by the solution curve $\mathcal Z$ defined in Lemma 3. Observe that even though the biologically relevant domain for (2.1) is just $\Omega$, $\mathcal{T}_{{\mathbb V}(\bar I, \bar u)}(x)$ may take finite values for points $x\notin\Omega$, while $\mathcal{T}_{{\mathbb V}(\bar I, \bar u)}(x) = +\infty$ for $x:x_1, x_2 < 0$.

For the purpose of our study we are interested in minimal entry times not exceeding a given time horizon $T$. As $\Omega$ is positively invariant under $f$, we employ the penalization term

To achieve higher accuracy of approximation for $w(t, x), t\in[0, T]$ close to the boundary $\partial\Omega$ we solve numerically (4.5) on a uniform stencil for a domain which contains $\Omega$ in its interior using a second-order Runge-Kutta scheme for the time integration. We find an approximation of the minimal entry time function on $\Omega$ by using the obtained solution values for $w(t, x), x\in\Omega$, as well as the optimal control using the reconstruction algorithm given in references ^[31,32].

5.   Results and discussion

Mathematical analysis of models for vector-borne diseases tends to concern the asymptotic rather than the transient behavior of their solutions using the basic reproduction number ${\mathcal R}_0$ as a criterion for the stability of the disease-free equilibrium ^{[34,35,36,37]}. Of interest are the parameter values that keep ${\mathcal R}_0$ below unity. Applications of optimal control theory often focus on optimal allocation of resources (by minimizing cost functionals of different form) without addressing future transient dynamic behavior ^[38].

In this work we study the latter question by considering a control problem with restricted state space, defined by a cap on the size of the infected host compartment, and address two major questions. In Question Q.1 we seek the viability kernel ${\mathbb V}$ for the model (2.1), the subset of the state space $\Omega$ that comprises those initial data such that a viable control $u(t)$ exists to maintain the size of infected host compartment $x_1(t)$ below the infection cap $\bar I$ for all future times. Then, we address a reachability problem formulated as Question Q.2. If the interventions start from an initial condition $x^0\notin{\mathbb V}$ which is outside the viability kernel, the objective stated in Question Q.1 cannot be met, so instead we are interested in finding the optimal control strategy $u(t)$ which reduces $x_1(t)$ below $\bar I$ in the shortest time. The motivation behind Question Q.2 relates to the fact that control measures studied here are restricted in duration due to technological and design limitations leading to a loss of product's repellent property.

Analysis of viability and reachability allows decision makers to assess the effect of this control intervention on future transient behavior given the epidemic situation at the start. This framework gives an opportunity to analyze the dynamic behavior of the controlled model on its entire domain of definition, or large subsets thereof, by describing existence of viable controls and of sets reachable in finite time. In particular, it is possible to estimate whether under the constraints of cost $C_{max}$ (which determines the population coverage $\bar u$) and duration of repellent action $T$, the objective of reducing the size of the infected host compartment below the infection cap $\bar I$ can be met.

The reachability analysis could serve as a predictive indicator for the performance of the intervention with parameters $\bar u, \bar I, T$ by providing analysis of the model's trajectories. If a trajectory which reaches the target (viability kernel) within the time interval $(0, T)$ exists, the campaign is considered effective. Else, the intervention strategy is to be considered insufficient and may need to be supplemented by additional measures (use of repellents with higher efficacy, use of LLINs, indoor residual spraying, vaccinations, preventive treatment, etc.).

We use parameter values for a malaria model fitted for Botswana and Zimbabwe from reference ^[19] (listed in Table 1) as a basis for the numerical simulations to exemplify this method. Viability kernels and minimal entry time functions are computed and compared across different values of $\bar u$.

The main differences between the parameters is the mosquito biting rate, $a_m = 0.082$ for Botswana and $a_m = 0.241$ for Zimbabwe, and the expected life-time of a mosquito $1/\mu$, which is three times higher in Botswana. According to reference ^[19] and therein cited references, the parameter $\mu$ reflects a difference in overall coverage of indoor residual spraying in the two countries. This intervention influences the vector's demography and ecology by direct insecticidal or repellent action, but different spraying solutions are used in homes built from traditional materials (wood, clay or mud bricks) and in Western-style housing ^[39]. While it has suppressed the transmission of malaria in the past, societal refusal to indoor spraying due to pungency, home wall staining, bug infestation, etc. is on the rise ^[40], and the positive trends in disease control have been reversed.

The results from Proposition 1 enable us to compute the threshold for population coverage $\min\bar u$ such that for a given infection cap $\bar I$ the viability kernel ${\mathbb V}(\bar I, \bar u)$ is not trivial:

If $\min\bar u\le0$, then the viability kernel ${\mathbb V}(\bar I, \bar u)$ has a positive Lebesgue measure in $\mathbb{R}^2$ for all $0\le\bar u\le1$.

Figure 1 plots the trade-off relationship (5.1) between $\bar I, \min\bar u$ for both countries for two values of repellent efficacy $k$. It shows the minimal coverage $\bar u$ required to produce a viability kernel of positive Lebesgue measure for the given infection cap $\bar I$. An immediate consequence is that in the case of Botswana, for $\bar I\ge 0.01$, any extent of the control measures will suffice to produce a non-trivial viability kernel. In contrast, more control efforts must be employed in the case of Zimbabwe, where at least 20% population coverage is required ($\min \bar u > 0.2$) to produce a non-trivial viability kernel. Of course, we must look at the area of the viability kernel to evaluate better the potential for keeping the epidemic peak below $\bar I$ with the repellent-based strategy under consideration.

The infection cap is $\bar I = 0.1$ and repellent efficacy $k = 0.6$. The epidemiological parameters for Botswana indicate that the viability kernels are not trivial for all $0\le\bar u\le 1$. Figure 2a) displays the area of the viability kernel for a range of values of $\bar u$. There is a positive correlation between $\bar u$ and the area of ${\mathbb V}(\bar I, \bar u)$. For the given epidemiological parameters in Zimbabwe, the situation is different: for $\bar I = 0.1$ the viability kernel is trivial if $\bar u < 0.2152$; in other words, under this assumption the model predicts it is not possible to maintain the size of infected host compartment below $\bar I$ for any initial data in $\Omega$ if the population coverage is less than 21.5%. The area of the viability kernel for a range of values of $\bar u\ge 0.215$ is plotted in Figure 2b).

Examples of viability kernels' shape for Botswana are shown in Figure 3. The viability kernels for Zimbabwe have a similar shape as those for Botswana, displayed in the panels of Figure 3, and are not shown.

An example for a reconstructed optimal trajectory inside $\Omega$ is given in Figure 4a), and the control effort – in Figure 4b) – it is computed using the algorithm from reference ^[31]. Note that for the same initial data, without intervention, the size of infected host compartment $x_1$ exceeds the cap $\bar I$ (the trajectory follows the purple curve in Figure 4a).

The numerical approximations for the minimal entry time functions $\mathcal{T}_{{\mathbb V}(\bar I, \bar u)}$ are shown in Figure 5 for Botswana and Figure 6 for Zimbabwe. In the simulations we use an infection cap of 10% of the population ($\bar I = 0.1$) and a time horizon of three months ($T = 100$ days). The results show that with these constraints, under the parameter values for Botswana, an outbreak would be successfully brought under the infection cap within the given time horizon, unless the proportion of infected mosquitos is very high. In fact, for $\bar u = 0.6$, the set of initial data $x^0$ such that ${\mathcal T}_{{\mathbb V}}(x^0) < 100$ days (those such that the size of the infected host compartment can be reduced to less than 0.1 within 100 days) is contained within the rectangle $[0, 1]\times[0, 0.25]$ and for $\bar u = 0.9$, it coincides almost fully with the rectangle $[0, 1]\times[0, 0.5]$. That means that increasing the coverage from 60% to 90% of the population approximately doubles the area of the set $\{x^0\; |\; {\mathcal T}_{{\mathbb V}}(x^0) < 100\mbox{ days}\}$. However, at very high levels of infected mosquitos at the start of the intervention ($x_2(0) > 0.5$) the intervention cannot be considered effective as for all different extents of coverage $\bar u\in[0.6, 0.9]$ the minimal time needed to reduce the outbreak size below the prescribed cap exceeds 100 days.

For Zimbabwe, the minimal entry time functions show that for the same maximum population coverage with repellent $\bar u$ less time is required to bring the outbreak under the chosen infection cap $\bar I$. The isolines in Figure 6 follow a radial pattern rather than parallel shifts like those for Botswana. Even at $\bar u = 0.6$ for the entire set of initial conditions $x^0\in[0, 1]\times[0, 0.5]$ the minimal entry time to the viability kernel is less than 90 days. Reconstructed optimal trajectories for the reachability problem in Question Q.2 for an initial condition $x^0 = (0.3, 0.2)$ are plotted in Figure 7.

It may seem paradoxical that for the model with parameters for Zimbabwe the minimal entry time to the respective viability kernel is shorter than for Botswana; in other words, a scenario with a higher basic reproduction number ${\mathcal R}_0$ appears more elastic in its response to the campaign of repellent-based interventions modeled by (2.1). That elasticity could be explained with the difference in the vector's ecology: the model is parameterized with a threefold higher mortality rate for the mosquitos in Zimbabwe ($\mu = \frac1{10}$) relative to Botswana ($\mu = \frac1{30}$), despite the fact that the pathogen transmission rate from vector to host $\alpha$, and ${\mathcal R}_0$ for a baseline scenario of no control ($u(t) = 0$) are higher for Zimbabwe ^[19]. This simple example shows how transient dynamical properties of the controlled system (2.1) in the case of endemicity stay decoupled from the value of a traditional epidemiological metric ${\mathcal R}_0$.

These simulations agree with {entomological} observations of the importance of interventions which target the vector's demography such as indoor residual spraying, LLINs, etc, for suppressing pathogen transmission ^[41]. Different societal behavior and inconsistent following of measures such as indoor spraying ^[19] can potentially have an impact on the performance of interventions such as those considered here and their ability to suppress sporadic outbreaks in the shortest time. The extent of compliance with a measure that is considered primary for malaria control by the WHO may influence the performance of alternative interventions reliant on {repellent-based personal protection}. This result stresses the importance for a multi-modal approach in controlling vector-borne disease outbreaks.

Our model has several limitations: first, it assumes that the host population is homogeneous, and does not consider age structure or differentiation in terms of disease status (such as asymptomatic, mild, severe and in need of treatment). Second, the host and vector populations are taken constant in time, ignoring effects of seasonal changes on the transmission dynamics. Third, it assumes that all hosts receive the same level of protection by the control effort, but population heterogeneity due to variability in repellent action may exist. Finally, (2.1) does not model immunity in the population. These issues can be addressed by increasing model complexity, including seasonal variations in the mosquito population, adding more compartments for the host population, etc. The state constraint (infection cap) can also be redefined to the respective class of infected (e.g., those in need of treatment). The proposed analytical framework can be adapted for extensions of the model (2.1), such as a time-delay system, a combination of measures (use of insecticides or vaccinations), or a system with an infected and infectious vector compartments ^[42].

6.   Conclusion

We have used optimal control theory to solve problems of existence of viable controls and of reachability for a simple model of a vector-borne disease. Available mathematical tools exist that allow decision and policy makers to assess effects of the control intervention on transient dynamics of the model in finite times, rather than just to determine asymptotic properties of equilibria. This approach reveals that, sometimes, depending on the epidemiological parameters and the state of the epidemic (determined by the size of infected host and vector compartments), even extending control efforts to the maximum may not suffice to meet policy objectives such as reduction of the size of the infected compartment below a given level in a fixed time span. Hence, existence of controls that satisfy the constraints of the intervention can be interpreted as an efficiency metric, which suggest whether the current intervention is sufficient or a combination of additional intervention modes may be required. Of course, to be able to make reliable predictions for the performance of a given control intervention, the model should use parameters rightful for the provided epidemiological and/or entomological data.

The presented model is based on assumptions of population homogeneity, simple disease dynamics and linear cost function for the control, and bringing the model closer to real-world scenarios will undoubtedly increase its complexity. Adding additional state variables means adding an extra dimension to the Hamilton-Jacobi-Bellman equation, and increases the computational complexity for the numerical approximation of the value function and the minimal entry time function. This obstacle should be addressed by development of appropriate, more efficient numerical methods for solving the associated optimal control problem.

Acknowledgments

The author would like to thank Bob W. Kooi for helpful discussion during the preparation of the manuscript. This work is supported by the Bulgarian National Science Fund (FNI) within the National Scientific Program Petar Beron i NIE of the Bulgarian Ministry of Education [contract KP-06-DB-5]

Conflict of interest

The author declares there is no conflict of interest.

Supplementary Material

Proof of Lemma 2. For part (i), consider the Lyapunov candidate function on $\Omega$:

Note that on $\Omega\in{\mathbb R}^2_+$, ${\mathcal L}\ge0$ with equality attained only at the disease-free equilibrium at the origin $\mathcal{O}$.

Note that

everywhere on ${\mathbb R}^2_+\setminus\mathcal{O}$ because

Hence ${\mathcal L}$ is positive definite, $|{\mathcal L}|\to+\infty$ on ${\mathbb R}^2_+$ and $d\mathcal{L}/dt$ is negative definite on ${\mathbb R}^2_+\setminus\mathcal{O}$. Then Lyapunov's stability theorem [43,Chapter 5.26] implies that the disease-free equilibrium $\mathcal{O}$ is globally asymptotically stable.

For part (ii), the global asymptotic stability of the endemic equilibrium (3.5) for the 2-dimensional system (2.1) is shown using the Lyapunov candidate function:

Note that ${\mathcal L} > 0$ and $|{\mathcal L}|\to+\infty$ on $\Omega$ (for $y_1\to0, y_2\to0$). After some algebraic transformations we arrive to

As $y_1, y_2 > 0$, equality $d\mathcal L/dt = 0$ can be achieved only at $\mathcal{E}_\ast$. Thus ${\mathcal L}$ is positive definite and $d\mathcal L/dt$ is negative definite on $\Omega\setminus\mathcal{O}$. Lyapunov's stability theorem implies that the endemic equilibrium $\mathcal{E}_\ast$ is globally asymptotically stable.

Proof of Proposition 1. The proof of case (i) is based on the local asymptotic stability of the disease-free or the endemic equilibrium of the system (2.1) with constant controls. In particular, we distinguish between 2 cases: first, if $I_{crit}\le0$, then $\mathcal{O}$ is globally asymptotically stable for (2.1) with $u(t)\equiv \bar u$ (see Lemma 2). This means that trajectories starting from any neighbourhood of $\mathcal{O}$ that is sufficiently small and contained inside $\mathcal{A}(\bar I)$ will converge to the equilibrium, maintaining the constraint on $x_1(t) < \bar I$ for all $t > 0$.

Second, if $I_{crit} > 0$, then the endemic equilibrium $\mathcal{E}_\ast$ is globally asymptotically stable for (2.1) with $u(t)\equiv\bar u$, and $\mathcal{E}_\ast\in\mathcal{A}(\bar I)$ (see Lemma 2). To show that $\mathtt{int}\; \mathbb{V}$ has a positive Lebesgue measure in $\mathbb{R}^2$, we consider the lens-shaped set $\mathcal{D}_0$ bounded by the $x_1$- and $x_2$-nullcline inside $\mathcal{A}(\bar I)$, and show it is a viability domain for (2.1) (see Figure S1). The $x_1$-nullcline is given by

with the outward-pointing normal

For $x\in\mathcal{N}_1$, it holds $f_2(x, \bar u) = \frac{x_1(\alpha\beta\nu(1-k\bar u)^2-\gamma\mu)- x_1^2(\alpha\beta\nu(1-k\bar u)^2+\gamma\beta(1-k\bar u))}{\alpha(1-k\bar u)\nu(1-x_1)}\ge 0$ with equality at $\mathcal{O}$ and $\mathcal{E}_\ast$, so $\langle f(x, \bar u), \bf{n}_{\mathcal{N}_1}\rangle = -\frac{\gamma f_2(x, \bar u)}{\alpha(1-k\bar u)\nu(1-x_1)^2}\le0$, and the curve (6.1) is impermeable from $\mathtt{int}\; \mathcal{D}_0$ for the flow $f(x, \bar u)$. The $x_2$-nullcline is given by

with the outward-pointing normal

For $x\in\mathcal{N}_2$, it holds $f_1(x, \bar u) = \frac{x_1(\alpha\beta\nu(1-k\bar u)^2-\gamma\mu)- x_1^2(\alpha\beta\nu(1-k\bar u)^2+\gamma\beta(1-k\bar u))}{\beta(1-k\bar u)x_1+\mu}\ge0$ again with equality at $\mathcal{O}$ and $\mathcal{E}_\ast$, so $\langle f(x, \bar u), \bf{n}_{\mathcal{N}_2}\rangle = - f_1(x, \bar u)\le0$ and the curve (6.2) is impermeable from $\mathtt{int}\; \mathcal{D}_0$ for the flow $f(x, \bar u)$. We conclude that this lens-shaped set $\mathcal{D}_0$ is forward invariant under $f(x, \bar u)$ and $\mathcal{D}_0\subset\mathcal{A}(\bar I)$ is a viability domain, completing the proof of case (i).

The proof of (ii) uses the componentwise inequality $f(x, u) \ge f(x, \bar u), \forall x\in\Omega$ to invoke the comparison theorem ^[20] and obtain

so $\bf{x}^{\bar u}(t; x^0)$ is a subsolution for (2.1) for all $x^0$. Lemma 2 establishes that for any initial condition $x^0\in\Omega\setminus\mathcal{O}$, this subsolution converges to $\mathcal{E}_\ast\in\Omega\setminus\mathcal{A}(\bar I)$. In other words, there exists $\tau > 0$ such that for the subsolution $\bf{x}^{\bar u}(t; x^0)$, $x_1^{\bar u}(t) > \bar I$ for all $t > \tau$. The comparison (6.3) implies that for the any feasible trajectory $\bf{x}^u(t; x^0)$ of the control system (2.1), $x_1^{ u}(t)\ge x_1^{\bar u}(t) > \bar I$ for all $t > \tau$. Thus, the only forward invariant set inside $\mathbb{V}(\bar I, \bar u)$ is the disease-free equilibrium $\mathcal{O}$.

Proof of Lemma 3. Denote for the sake of shortness $\tilde \alpha = \alpha\nu(1-k\bar u), \tilde \beta = \beta(1-k\bar u)$, and $\tilde f_1(z, s) = f_1(z, s, \bar u), \tilde f_2(z, s) = f_2(z, s, \bar u)$, and $H(s, z) = \frac{\tilde f_1(z, s)}{\tilde f_2(z, s)}$ the right-hand side of (3.6).

Note that due to the choice of $\bar v$, the numerator $\tilde f_1$ of $H$ is 0 at $s = \bar v$. Further, the denominator of $H$ at $(\bar v, \bar I)$ satisfies

so long as $\tilde\alpha\tilde\beta-\gamma\mu\le0$ (when $I_{crit}\le0$) or when $\bar I > I_{crit}, \tilde\alpha\tilde\beta-\gamma\mu > 0$. Therefore,

We have

and in particular, $\tilde f_1(z, s) > \tilde f_1(\bar I, \bar v)$ and $\tilde f_2(z, s) < \tilde f_2(\bar I, \bar v)$ if $s > \bar v, z < \bar I$. Also we have $\max \tilde f_1(z, s) < +\infty, \min \tilde f_2(z, s) > -\infty$ on $\Omega$.

This implies that the denominator of $H$ is negative and bounded away from 0 on $[\bar v, 1]\times [0, \bar I]$ due to (6.4). Consequently, being a rational function, $H$ is continuous in $s$ and uniformly Lipschitz in $z$ on $[\bar v, 1]\times [0, \bar I]$. By the Picard-Lindelöf theorem, there exists $\varepsilon > 0$ such that (3.6) has a unique $C^1$-solution $z(s)$ for $s\in[\bar v, \bar v+\varepsilon)$. We show now that this solution $z(s)$ exists on the entire interval $[\bar v, +\infty)$.

We verify that $z(s)$ is monotone decreasing on $[\bar v, \bar v+\varepsilon)$ based on (6.5). Since $z$ is $C^1$, its Taylor expansion around $\bar I = z(\bar v)$ reads

Therefore, using the relation $\tilde f_1(\bar I, \bar v) = 0$:

Thus, there exists $\varepsilon > 0$ such that for values $\bar v < s < \bar v+\varepsilon$, $\tilde f_1(z(s), s) > 0$, while $\tilde f_2(z(s), s) < 0$ due to (6.4). Hence, the solution $z$ is monotone decreasing on $[\bar v, \bar v+\varepsilon)$, and due to the monotonicity of $\tilde f_1, \tilde f_2$ provided by (6.6), $H < 0$ strictly on $(\bar v, +\infty)\times[0, \bar I)$. Because the solution $z(s), s > \bar v$ is bounded by $z(\bar v) = \bar I$ irrespectively of $\varepsilon$, it can be continued on the whole interval $[\bar v, +\infty$, and remains monotone decreasing because $dz/ds < 0$.

Note that it might happen that $z(v^\ast) = 0$ for some $\bar v < v^\ast\le1$. Then we restrict the solution curve $\mathcal{Z}$ to the domain $[\bar v, v^\ast]$. This completes the proof.

Proof of Proposition 2. The proof follows the ideas from reference ^[17]. To show that the proposed sets are viability kernels, we have to verify their forward invariance under $f(x, u)$ (that is, that they are viability domains for $f$) and their maximality.

Case (i). Consider the segment $\bf{s} = \{\bar I\}\times[0, 1]$, which forms the eastern boundary of $\mathcal{A}(\bar I)$ in the $(x_1, x_2)$-state space and the outward-pointing normal to $\bf s$, $\bf{n}_s = (1, 0)^T$. Then $\langle f(\bar I, x_2, \bar u), \bf{n}_s\rangle = \alpha (1-k\bar u)\nu x_2 (1-\bar I)-\gamma \bar I\le(1-k\bar u) \alpha\nu (1-\bar I) -\gamma \bar I < 0$. This means there exists a control function $u\in\mathcal{U}$ such that the entire set $\mathcal{A}(\bar I)$ is forward invariant under (2.1), and therefore, $\mathbb{V}(\bar I, \bar u) = \mathcal{A}(\bar I)$.

Case (ii). Let $\mathcal{D} = \{[0, \bar I]\times[0, \bar v]\}\cup \{ (x_1, x_2) | \bar v\le x_2\le\min\{v^\ast, 1\}, 0\le x_1\le z(x_2)\}$. Recall $\bar v$ from Lemma 3. Observe that on the segment $\bf s_2 = \{\bar I\}\times(\bar v, 1]$, the inner product of the velocity field $f$ and the outward-pointing normal $\bf n_2$ to $\bf s_2$, $\langle f(\bar I, x_2, u), \bf{n}_2\rangle = \alpha (1-ku)\nu x_2 (1-\bar I)-\gamma \bar I > (\frac{1-ku}{1-k\bar u}-1)\gamma \bar I\ge0$. Hence, by continuity of $f$, a neighbourhood of $\bf s_2$ in $\mathcal{A}(\bar I)$ is not contained inside $\mathbb{V}(\bar I, \bar u)$. On the other hand, on the segment $\bf s_1 = \{\bar I\}\times[0, \bar v]$ the inner product of the velocity field $f$ and the outward-pointing normal $\bf n_1$ to $\bf s_1$, $\langle f(\bar I, x_2, \bar u), \bf{n}_1\rangle\le 0$. Finally, the outward-pointing normal to $Z$ is $\bf{\tilde n} = (1, -z'(x))^T$. Then $\langle f(z(x_2), x_2, \bar u), \bf{\tilde n}) = 0$. Hence, by choosing $u = \bar u$, $\mathcal Z$ is impermeable to the flow of $f(x, \bar u)$. Thus, there exists $u$ such that the interior $\mathtt{int}\mathcal{D}$ is forward invariant under $f(x, u)$.

Finally to conclude that $\mathcal{D}$ is forward invariant, we show that the curve $\mathcal Z\in\mathcal D$, itself being a feasible trajectory for (2.1) with $u = \bar u$. Indeed, the inequality $x_2^0 > \bar v, x_1^0 < \bar I$ holds for any initial data $x^0 = (x_1^0, x_2^0)\in \mathcal{Z}$. Therefore, $\tilde f_2(x_1^0, x_2^0) < \tilde f_2(\bar I, \bar v) < 0$ by (6.4), and since $\frac{d}{dt}x_2(x^0) < 0$, it follows $x_2(t) < x_2^0, t > 0$.

Let $\tau_{\mathcal{B}}(x^0)$ be the minimal entry time function (4.1) to the target set $\mathcal{B} = \{(x_1, x_2) | \ x_1 > \bar I \vee x_2 < \bar v \}$ for the control $u(t) = \bar u$, or

Since $\mathcal{B}$ contains the globally asymptotically stable endemic equilibrium $\mathcal{E}_\ast$ (3.5), $\tau(x^0) < +\infty$ for all $x^0\in \mathcal Z$.

We choose any initial condition $x^0\in\mathcal Z$. Then the trajectory of (2.1) for $u = \bar u$ satisfies $x_2(t) < \bar v$ for $t > \tau(x^0)$. Thus, $\{x_2(t)|t\in[0, \tau(x^0)]\}\subset[\bar v, \min\{ v^\ast, 1\}]$ and $x_2(t)$ takes values inside the domain of the function $z$ which solves the problem (6.4). On $t\in[0, \tau(x^0)]$ the difference between $x_1(t)$ and $z(x_2(t)$ remains constant in $t$ because

Due to the initial condition $x^0\in \mathcal Z$, $x_1(0)-z(x_2(0)) = 0$, so the above implies $x_1(t)-z(x_2(t)) = 0, t\in(0, \tau(x^0))$ implying that the curve $\mathcal Z$ is part of a feasible trajectory for (2.1).

As it holds by construction, either $x_1(\tau(x_0)) = \bar I$ or $x_2(\tau(x_0)) = \bar v$, by the inverse function theorem for $z$ we obtain that $\bf{x}_{x^0}^{\bar u}(\tau(x^0)) = (\bar I, \bar v)$. Since $x^0\in\mathcal Z$ was chosen arbitrary, the entire curve $\mathcal Z$ as defined by the Lemma consists of a trajectory for (2.1). We have demonstrated that the set $\mathcal D$ as defined in (3.6) is closed, forward invariant under $f(x, u)$ and it is a viability domain for $f(x, u)$.

It remains to show that $\mathcal D$ is the maximal viability domain for $f(x, u)$. Consider an initial condition $y^0\in\mathcal{A}(\bar I)\setminus \mathcal D$ for $f$ and {claim} the trajectory $\bf{x}^{\bar u}(t, y^0)$ leaves $\mathcal{A}(\bar I)$ in finite time. It holds $y_1^0 < \bar I$, but $y_2^0 > z^{-1}(y_1^0)$ due to the definition of $\mathcal D$. Denote $\tilde y^0 = (y_1^0, z^{-1}(y_1^0))$, meaning $\tilde y^0\le y^0$ componentwise. Since the map $f$ is quasimonotone, we can use the comparison principle [24,Chapter 3] to obtain the componentwise inequalities

Therefore, the first element in the difference vector $\Delta(t) = \bf{x}^{\bar u}(t; {y^0})- \bf{x}^{\bar u}(t; {\tilde y^0})$, which we denote by $\Delta_1(t)$ satisfies $\Delta_1\in C^1(0, +\infty), \Delta_1(t)\ge 0$.

Denote $\tau_0 = \tau(\tilde y^0)$, which means $\tilde y_1(\tau_0) = [\bf{x}^{\bar u}(\tau_0;\tilde y^0)]_1 = \bar I$. Assume $y_1(\tau_0) = [\bf{x}^{\bar u}(\tau_0;y^0)]_1 = \bar I$, which translates to $\Delta_1(\tau_0) = 0$, and $\Delta_1$ having a local minimum at $t = \tau_0$. Hence,

Plugging in the equality $y_1(\tau_0) = \tilde y_1(\tau_0) = \bar I$ into (6.7), we see that

and obtain $\bf{x}^{\bar u}(\tau_0;y^0) = \bf{x}^{\bar u}(\tau_0;\tilde y^0)$. With $y\neq \tilde y$, this equality is a contradiction to the uniqueness of the solution trajectory to the autonomous system $\frac{d}{dt} y = f(y, \bar u)$. Therefore, $y_1(\tau_0) > \tilde y_1(\tau_0) = \bar I$, showing that $y^0\notin\mathbb{V}(\bar I, \bar u)$. This argument shows $\mathcal D$ is maximal, and we conclude $\mathbb{V}(\bar I, \bar u) = \mathcal D$.