A minimal model coupling communicable and non-communicable diseases

M. Marvá∗, E. Venturino∗∗, M.C. Vera∗

∗ Universidad de Alcalá, Departamento de Física y Matemáticas.

Member of the research group Nonlinear Dynamics and Complex Systems.

Alcalá de Henares, Spain.

∗∗Dipartimento di Matematica ”Giuseppe Peano”.

Member of the INdAM research group GNCS.

Univertitá di Torino, Torino, Italy.

###### Abstract

This work presents a model combining the simplest communicable and non-
communicable disease models. The latter is, by far, the leadingn cause of
sickness and death in the World, and introduces basal heterogeneity in
populations where communicable diseases evolve. The model can be interpreted
as a risk-structured model, another way of accounting for population
heterogeneity.

Our results show that considering the non-communicable disease (in the end,
heterogeneous populations) allows the communicable disease to become endemic
even if the basic reproduction number is less than $1$. This feature is known
as subcritical bifurcation. Furthermore, ignoring the non-communicable disease
dynamics results in overestimating the reproduction number and, thus, giving
wrong information about the actual number of infected individuals. We
calculate sensitivity indices and derive interesting epidemic-control
information.

Keywords: Non-communicable disease, communicable disease, basic reproduction
number, subcritical bifurcation, supercritical bifurcation, heterogeneous
populations, risk-structured. PACS: 34, 37N25, 92D30

##  1 Introduction

It is a fact that non-communicable diseases (NCDs) include cardiovascular
disease, cancer, chronic respiratory disease, or diabetes, are the main cause
of sickness and death worldwide [26]. Just in 2000 NCDs were responsible for
35 million deaths (about 60% of all deaths around the world), and in 2020
these figures raised to 41 million death (i.e. about 71% of deaths worldwide)
[26]. NCDs are the result of a combination of non-reversible genetic and
physiological factors but also of environmental and behavioral factors that
may be reverted. Use of or exposure to tobacco, alcohol abuse, unhealthy
diets, or physical inactivity [7] are among these revertible factors, as well
as air pollution and environmental contamination [22]. NDCs are very common
and, therefore, play a key role in the epidemiology of communicable or
infectious diseases (CD) [9], [27].

The reproduction number, $R_{0}$, is a key quantity in the dynamics of a
communicable disease. $R_{0}$ stands for the average number of secondary
infections produced by an infected individual in a population made just of
susceptible individuals [2]. It is well known that $R_{0}>1$ enables
communicable diseases to become endemic. However $R_{0}<1$ does not always
lead to the eradication of the communicable disease. This somewhat counter-
intuitive fact is known as subcritical bifurcation [12] (also often less
properly named backward bifurcation [15]). Mechanisms leading to subcritical
bifurcations in epidemiological models are proposed in [13], while general
necessary and sufficient conditions for an epidemiological model to display a
subcritical bifurcation are obtained in [5]. This phenomenon has important
consequences from the viewpoint of epidemics control, since reducing $R_{0}$
below $1$ may not be sufficient to avoid the disease endemic scenario. An
instance of this unfortunate behavior is shown in the case of TB in India,
[25].

This work is aimed to analyze the interplay between NCDs and CDs. For this
purpose, in Section 2 we set up a minimal model using one of the simplest
transmission laws [2] for CDs and the minimal number of epidemiological
stages. The structure of the model presented herein can be seen as a
simplified variant of risk-structured SIS models [16], [13], [17]. In doing
so, we isolate the net effect of the NCD/risk-structure on the behavior of the
CD. Thus, we disentangle the role of simple heterogeneity in the screened
population from other processes (see the discussion in Section 4). In Section
3 we analyze the model and derive sufficient and necessary conditions enabling
a subcritical bifurcation. We discuss the results and its implications on the
control of the CD in Section 4.

##  2 Model formulation

We focus on the interplay between the CD and the NCD. Thus, we seek the
minimal settings avoiding, for instance, demography processes. This approach
yields a laboratory model that can be expanded to face more general settings.

We assume that individuals affected by the NCD are somehow weaker to face the
CD. At time $t$ the population is partitioned into susceptible individuals
$S(t)$, i.e. those that are affected neither by the CD nor the NCD, weakened
individuals $W(t)$, i.e. those that suffer from the NCD but are not infected
by the transmissible disease, and individuals infected by the CD, $I(t)$,
regardless of whether they are weakened or not.

The model is built with ordinary differential equations. Next, we define the
rates at which individuals move from one compartment to each other, see Figure
1. We assume that susceptible individuals become weakened at a constant rate
$a$ and get out of the weakened class also at a constant rate $b$. It may
refer to many different processes. For instance, for diseases not directly
transmissible through contact (for instance, smoking) still the “observation”
of the other people’s opposed behaviors may influence both susceptible to take
the same habit, or conversely to suggest “addicted” people (i.e. $W$
individuals) to leave stop smoking. On the other hand, external dynamics may
impoverish/take out of poverty people, which is also a risk factor in front of
a CD. Note that the transmission of the NCD could still be formulated
considering more sophisticated functional terms [23], [21], but we disregard
this approach here to keep the model formulation minimal as mentioned. The
dynamics of the transmissible disease is somehow similar to (but not exactly)
the classical SIS model [2] with density-dependent transmission. Susceptible
and weakened individuals are assumed to behave differently concerning the CD
so that the transmission rates $\beta_{S}$ and $\beta_{W}$ differ from each
other. We assume that infected individuals recover and become immune at rate
$\gamma$. Note that we do not know whether an infected individual suffers or
not from the NCD. Thus, we do not care about recovered individuals (that leave
the model). We further consider that susceptible (in front of the CD)
individuals are introduced at rates $\gamma_{S}I$ and $\gamma_{W}I$ in the
corresponding compartment $S$ and $W$, so that the population size is kept
constant.

The analysis can be done mutatis mutandi also by considering frequency-
dependent transmission.

Figure 1: Flow diagram associated to system (1).

The ordinary differential equations system produced by the above-stated
hypotheses consists of two coupled submodels: one describing the communicable
disease and another one that describes the non-communicable disease. The
combined model reads as follows

$\left\\{\begin{array}[]{l}S^{\prime}=-aS+bW-\beta_{S}SI+\gamma_{S}I,\\\ \\\
W^{\prime}=aS-bW-\beta_{W}WI+\gamma_{W}I,\\\ \\\
I^{\prime}=\beta_{S}SI+\beta_{W}WI-\gamma I,\\\ \end{array}\right.$ (1)

where $\gamma_{S}+\gamma_{W}=\gamma$. Note that the total population size
$N(t):=S(t)+W(t)+I(t)=N$ remains constant, as $N^{\prime}(t)=0$.

For sake of completeness, let us revisit these two well known models: the SIS
and the model for non-communicable diseases.

##### SIS model (CD).

If susceptible individuals are all of the class, system (1) simplifies into

$\left\\{\begin{array}[]{l}S^{\prime}=-\beta_{S}SI+\gamma_{S}I,\\\ \\\
I^{\prime}=\beta_{S}SI-\gamma_{S}I,\\\ \end{array}\right.$ (2)

where $S(t)+I(t)=N$ remains constant over time. It is nothing but the
classical SIS model with density-dependent transmission [2]. A straightforward
analysis reveals that system (2) possesses two equilibrium points: the trivial
equilibrium $(S^{*},I^{*})=(N,0)$ (no infected individuals) and a the endemic-
disease equilibrium

$(S^{*},I^{*})=\left(\dfrac{\gamma_{S}}{\beta_{S}},\,N-\dfrac{\gamma_{S}}{\beta_{S}}\right)$
(3)

The well known basic reproduction number

${R}_{0,S}=\dfrac{\beta_{S}}{\gamma_{S}}N$ (4)

determines whether the disease-free equilibrium (${R}_{0,S}<1$) or the
endemic-disease equilibrium (${R}_{0,S}>1$) is the global attractor of system
(2).

##### Non-communicable disease.

In the absence of the infectious disease ($I(t)=0$), system (1) reduces to
system

$\left\\{\begin{array}[]{l}S^{\prime}=-aS+bW,\\\ \\\ W^{\prime}=aS-bW,\\\
\end{array}\right.$ (5)

where $S(t)+W(t)=N$ remains constant over time. We will discuss in section 4
further extensions of this submodel. We assume $N>0$ (there are individuals)
so that there exists a non trivial equilibrium point

$(S^{*},W^{*})=\left(\dfrac{b}{a+b}N,\,\dfrac{a}{a+b}N\right)$ (6)

where $b/(a+b)$ and $a/(a+b)$ are the fraction of susceptible and weakened
individuals within the entire population $N$. It is straightforward that the
nontrivial equilibrium (6) is a global attractor. This feature can be
interpreted as the non-communicable disease being structural to the
population.

##  3 Results

In this section, we analyze the long term behavior of the solutions of system
(1), that is, the so-called equilibrium points and their stability.

A first step consists of showing that the model is well behaved, that is,

###### Proposition 1

The solutions of system (1) are bounded from above and the non negative cone
is forward invariant.

Proof: All the solutions of system (1) are bounded since the total population
size is kept constant because $S^{\prime}(t)+E^{\prime}(t)+I^{\prime}(t)=0$.
The invariance of the non negative cone

$\bar{\mathbb{R}}_{+}^{3}:=\left\\{(S,W,I)\in\mathbb{R}^{3};\,S\geq 0,\,E\geq
0,\,I\geq 0\right\\}$

is equivalent to prove that any solution with initial values on the boundary
can not become negative as the time flows. For instance, assume that
$W(t_{0})=0$ and $S(t_{0})\cdot I(t_{0})\neq 0$. It follows that
$W(t_{0})^{\prime}=aS(t_{0})+\gamma_{W}I(t_{0})>0$, so that the corresponding
solution grows towards positive values, and the same holds assuming
$S(t_{0})=0$, $W(t_{0})\cdot I(t_{0})\neq 0$.

Assume now that $I(t_{0})=0$ and $S(t_{0})\cdot W(t_{0})\neq 0$. Then
$I^{\prime}(t_{0})=0$ regardless of the value of $S(t_{0})$ and $W(t_{0})$.
There are no infected individuals and there will be none. The solution of the
system evolves constrained by $S(t_{0})+W(t_{0})=N=cte$, that is to say that
system (1) is reduced to system (5) and its solution converge to (6).
$\blacksquare$

As we have already said, the non-communicable disease is supposed to be
inherent to that population. The first result consists of determining
conditions so that an outbreak of the communicable disease lead to an endemic
disease scenario. In other words, we seek conditions enabling the semitrivial
equilibrium point

$E_{0}^{*}=(S_{0}^{*},W_{0}^{*},I_{0}^{*})=\left(\dfrac{b}{a+b}N,\,\dfrac{a}{a+b}N,\,0\right)$
(7)

to be asymptotically stable (communicable-disease-free state) or unstable
(endemic communicable disease scenario). Note that trivial equilibrium point
$E_{0}^{*}$ consists of the components of the nontrivial equilibrium (6) of
system (5) along with $0$ infected individuals in the third entry.

###### Proposition 2

Consider system (1), the disease free equilibrium $E_{0}^{*}$ given by (7) and
the reproductive number $R_{0}$ defined by

$R_{0}=\dfrac{b\beta_{S}+a\beta_{W}}{(a+b)(\gamma_{S}+\gamma_{W})}N.$ (8)

Then $E_{0}$ is unstable if $R_{0}<1$ and locally asymptotically stable if
$R_{0}>1$.

Proof: It follows from a standard analysis of the sign of the eigenvalues of
the Jacobian matrix of the flow of system (1). $\blacksquare$

Thus the CD free equilibrium (7) is unstable if $R_{0}>1$, which leads to a CD
endemic scenario. On the other hand, $R_{0}<1$ implies that the CD free
equilibrium (7) is locally asymptotically stable, meaning that any CD outbreak
will fade out (at least if the number of infected individuals is small
enough). These results provide us with valuable but incomplete information.
Namely:

  1. 1.

What are the conditions leading the CD free scenario to be globally
asymptotically stable? I.e., what are the conditions ensuring that CD free
scenario will be achieved regardless of the strength of a potential outbreak?

  2. 2.

What role does the NCD play in the dynamics of communicable disease?

In other words, we are interested in the structure of the set of the positive
equilibrium points of system (1) and its stability. Positive equilibrium
points are the component-wise positive solutions to system

$\left\\{\begin{array}[]{l}0=-aS+bW-\beta_{S}SI+\gamma_{S}I,\\\ \\\ 0=aS-
bW-\beta_{W}WI+\gamma_{W}I,\\\ \\\ 0=\beta_{S}SI+\beta_{W}WI-\gamma I,\\\
\end{array}\right.$ (9)

where $\gamma=\gamma_{S}+\gamma_{W}$. Let us assume that $I(t_{0})\neq 0$,
since otherwise system (9) reduces to system (5). Weakened individuals behave
differently from susceptible individuals in front of the CD, being plausible
$\beta_{W}\geq\beta_{S}$ and $\gamma_{W}\leq\gamma_{S}$ with at least one of
the inequalities being strict. We assume $\beta_{W}>\beta_{S}$ through the
manuscript, since $\gamma_{S}$ and $\gamma_{W}$ are not proper recovery rates.

Using the fact that the total population $S(t)+I(t)+W(t)=N$ is constant,
direct calculations yield that the number of infected individuals is the
solution of a quadratic polynomial equation

$\Psi(I)=\alpha_{2}I^{2}+\alpha_{1}I+\alpha_{0}=0$ (10)

where

$\begin{array}[]{c}\alpha_{2}=-\dfrac{\beta_{S}\beta_{W}}{\beta_{W}-\beta_{S}}<0,\\\
\\\
\alpha_{1}=-\dfrac{\beta_{S}(b+\gamma_{W})+\beta_{W}(a+\gamma_{S})-\beta_{S}\beta_{W}N}{\beta_{W}-\beta_{S}},\\\
\\\
\alpha_{0}=-\dfrac{(\gamma_{S}+\gamma_{W})(a+b)-(b\beta_{S}+a\beta_{W})N}{\beta_{W}-\beta_{S}}.\end{array}$
(11)

Note that $\beta_{W}>\beta_{S}$, since weakened individuals are weaker in
front of the CD. It is immediate calculate the number of infected individuals
$I^{*}$ (if any) at equilibrium. Then, $S^{*}$ and $W^{*}$, the number of
susceptible and weakened individuals at equilibrium can be calculated from
$I^{*}$ according to

$\begin{array}[]{c}I_{\pm}^{*}=\dfrac{-\alpha_{1}\pm\sqrt{\alpha_{1}^{2}-4\alpha_{0}\alpha_{2}}}{2\alpha_{2}},\\\
S_{\pm}^{*}=\dfrac{\beta_{W}}{\beta_{W}-\beta_{S}}\left(N-\dfrac{\gamma_{S}+\gamma_{W}}{\beta_{W}}-I_{\pm}^{*}\right),\\\
\\\
W_{\pm}^{*}=\dfrac{\gamma_{S}+\gamma_{W}-\beta_{S}S_{\pm}^{*}}{\beta_{W}}.\end{array}$
(12)

Note that condition $\alpha_{1}^{2}-4\alpha_{0}\alpha_{2}<0$ yields no
infected individuals, so that we assume from now on that
$\alpha_{1}^{2}-4\alpha_{0}\alpha_{2}\geq 0$. Direct calculations yield

$\alpha_{0}=0\Leftrightarrow R_{0}=1,\quad\alpha_{0}<0\Leftrightarrow
R_{0}<1,\quad\alpha_{0}>0\Leftrightarrow R_{0}>1$

along with

$\alpha_{1}=0\Leftrightarrow\Delta=1,\quad\alpha_{1}>0\Leftrightarrow\Delta>1,\quad\alpha_{1}<0\Leftrightarrow\Delta>1,$

where

$\Delta=\left(\dfrac{a+\gamma_{S}}{\beta_{S}}+\dfrac{b+\gamma_{W}}{\beta_{W}}\right)\dfrac{1}{N}$
(13)

It is not difficult to classify all the possible qualitatively different
outcomes in terms of sign of $\alpha_{0}$ and $\alpha_{1}$ (given that
$\alpha_{2}<0$). That is to say, in terms of the value of $R_{0}$ and $\Delta$
relative to $1$. Figure 2 sketches the interesting cases.

Figure 2:  Each panel displays the parabola defined by equation (10) in the
$I$-$I^{\prime}$ plane for different combinations of $R_{0}$ and $\Delta$.
Solid dots are the feasible equilibrium points that are asymptotically stable
(black) and unstable (gray).

Panels in Figure 2 display $I^{\prime}=\Psi(I)$; that is to say, the intercept
with the horizontal axis is the amount of infected individuals at equilibrium
(since $I^{\prime}=0$). Note that the sign of $I^{\prime}(t)$ determines
whether the number of infected individuals $I(t)$ increases ($I^{\prime}>0$)
or decreases ($I^{\prime}<0$), which yields the stability of the equilibrium
points. For instance, consider the right bottom panel: let $I(t,I_{0})$ be the
solution to $I^{\prime}=\Psi(I)$ such that $I_{0}=I(t_{0},I_{0})$, for $I_{0}$
larger than the most right equilibrium point (say $I^{*}$). It is apparent
$I^{\prime}(t_{0},I_{0})<0$, so that $I(t,I_{0})$ is strictly decreasing for
any $t\in(t_{0},t_{1})$. Also, $I^{*}<I(t,I_{0})$ for any $t\in(t_{0},t_{1})$.
Let $t_{1}$ define the maximal interval where $I(t,I_{0})$ is decreasing.
$I(t,I_{0})$ is bounded from below, strictly decreasing and continuous. Then,
there exists

$\hat{I}=\lim_{t\to\infty}I(t,I_{0}).$ (14)

The uniqueness of solutions precludes $\hat{I}<I^{*}$. Also, $\hat{I}>I^{*}$
implies that $I^{\prime}(t,I_{0})<\xi<0$ for all $t\in(t_{0},t_{1})$, that is
a contradiction with the existence of limit (14). Then, $\hat{I}=I^{*}$, which
concludes the proof.

The stability of all the other equilibrium points in Figure 2 follows
reasoning as before.

Essentially, three scenarios are possible: the global CD free scenario (left
and central panels in the first row of Figure 2), the global endemic CD
scenario (right panel in the first row and left and central panels at the
second row of Figure 2), and a third intermediate one that predicts either
endemic disease or disease free scenarios depending on the initial amount of
infected individuals (right panel of the second row of Figure 2). These
features are better shown with a bifurcation diagram, being $R_{0}$ the
bifurcation parameter, see Figure 3. In this context, $\Delta$ is the so-
called direction of bifurcation, so that $\Delta<1$ leads to a subcritical (or
backward) bifurcation (right panel in Figure 3) and $\Delta>1$ yields a
supercritical (or forward) bifurcation (left panel in Figure 3).

Figure 3:  Bifurcation diagrams plotting the total number of infected
individuals at equilibrium for different values of $R_{0}$. The left panel
displays a supercritical bifurcation and the right panel displays a
subcritical bifurcation. In gray unstable equilibrium points and in black
asymptotically stable equilibrium points. Parameter values: ,
$\beta_{W}=12.1$; $\beta_{S}=2.6$; $\gamma_{S}=8.4$; $\gamma_{W}=33.3$,
$a\in[0.001,10]$

For $R_{0}<1$ the solution $I^{*}=0$ is locally asymptotically stable, i.e.
the presence of a small number of infected individuals is not enough to
trigger an epidemic disease state, and the infected population will fade away
regardless of the value of $\Delta$. However, due to management decisions or
natural causes, the values of the parameters involved in the expression of
$R_{0}$ may change and increase $R_{0}$ so that it crosses the threshold value
$R_{0}=1$. In such a case, the solution $I^{*}=0$ is destabilized, which means
that the communicable disease becomes endemic even if there is a little
initial amount of infected individuals.

The quantity $\Delta$ plays a key role when $I^{*}=0$ is stable ($R_{0}<1$):
namely if $\Delta>1$ then $I^{*}=0$ is globally asymptotically stable, which
means that any epidemic outbreak will fade away regardless of the initial
number of infected individuals. On the contrary, $\Delta<1$ implies that there
exists $R_{0}^{*}<1$ such that for each $R_{0}^{*}<R_{0}<1$ there exist a
threshold number of infected individuals given by $I_{-}^{*}$ (see equation
(12)) such that beyond it the disease becomes endemic and stabilize at
$I_{+}^{*}$ (see equation (12)). For $\Delta<1$ and $R_{0}<R_{0}^{*}$
$I^{*}=0$ becomes globally asymptotically stable.

The stability of all the other equilibrium points in Figure 2 follows
reasoning as before.

##  4 Discussion

We next discuss our finding on system (1). Its apparent simplicity allows,
instead, to emerge relevant features. First of all, we account for the
possibility of subcritical bifurcation under minimal settings, see Section
4.1. We also analyze the consequences of not considering explicitly the NCD in
Section 4.2. Then, in Section 4.3, we focus on the effect of control
strategies (via modifying the coefficients of the system). In particular, we
calculate the sensitivity indices, we derive bounds for these indices, we
focus on the effect of modifying several coefficients at once, and we reveal
possible unexpected consequences when trying to handle disease outbreaks.

###  4.1 Subcritical bifurcation

Regardless of the approach (CD vs NCD or CD risk structure) system (1) may
undergo a subcritical bifurcation, which is against the $R_{0}$-dogma [24]
that states that $R_{0}>1$ leads to disease endemicity while $R_{0}<1$ induces
disease eradication. Usual causes of subcritical bifurcation are the use of
imperfect vaccine [4], [14], structured immunity [24] or exogenous re-
infection in TB disease [10].

In [13] several other biological or epidemiological mechanisms are proposed
such as vaccine-induced immunity waning at a slower rate than natural
immunity, disease-induced mortality in vector-borne diseases, and differential
susceptibility in risk-structured models (related to the latter, see [16] and
[17]).

Subcritical bifurcations can be found also in co-infection by an opportunistic
disease model [19], where two communicable diseases were considered (one of
them with saturating treatment rate [18]).

A salient feature of the model proposed and analyzed here lies in its ability
to undergo a subcritical bifurcation while not incorporating any of the above-
mentioned mechanisms. Thus, we show that subcritical bifurcations in
epidemiology are not such a rare occurrence. On the contrary, plain
heterogeneity in the CD susceptible class is enough to make a subcritical
bifurcation possible.

Thus, system (1) shows that subcritical bifurcations may occur in epidemic
models simply by considering the dynamics associated with a heterogeneous
population, which can be seen as a common factor in the above-mentioned
models.

###  4.2 What if the non-communicable disease is not explicitly considered?

Thus, let us assume that there is no weakened individuals compartment so that
the CD follows the simplest model (2). Even if we do not consider an explicit
compartment, the NCD is present in the population, and the coefficients of
system (2) must reflect in some way this fact. It is reasonable assuming that
the NCD-induced population heterogeneity will be captured by any reasonable
sampling procedure performed to estimate the coefficients of the model by
weighting the corresponding transmission ($\beta_{S}$ and $\beta_{W}$) and
recovery ($\gamma_{S}$ and $\gamma_{W}$) coefficients. This hypothesis is
equivalent to assuming that the dynamics associated with the NCD has already
achieved an equilibrium, which is the usual assumption when dealing with time-
scale systems [1] (also known as quasi-steady-state approximation [11]) that
has been used in co-infection by an opportunistic disease models [20], [19],
where both diseases are transmissible. Thus a fraction $b/(a+b)$ of the total
population is free of the NCD and the remaining fraction $a/(a+b)$ is not.
Then, to obtain a fair comparison, transmission and the recovery rates are set
to

$\beta_{S}\dfrac{b}{a+b}+\beta_{W}\dfrac{a}{a+b},\qquad\gamma_{S}\dfrac{b}{a+b}+\gamma_{W}\dfrac{a}{a+b},$
(15)

respectively, which yield the corresponding basic reproductive number:

$\widehat{R}_{0}=\frac{b\beta_{S}+a\beta_{W}}{b\gamma_{S}+a\gamma_{W}}N,\qquad\qquad(\widehat{S}^{*},\,\widehat{I}^{*})=N\left(\frac{1}{\widehat{R}_{0}},\,N-\frac{1}{\widehat{R}_{0}}\right)$
(16)

Direct calculations yield

$\frac{\widehat{R}_{0}}{R_{0}}=1+\frac{a\gamma_{S}+b\gamma_{W}}{b\gamma_{S}+a\gamma_{W}}.$
(17)

That is, the ratio (17) is always larger than $1$, implying that explicit
consideration of the NCD dynamics in the model does matter. Neglecting its
effect leads to overestimating the basic reproductive number and, thus, i)
thinking of an endemic disease scenario that may be not real and ii)
overestimating the number of infected individuals at equilibrium.

Figure 4 displays the bifurcation diagram of the total amount of infected
individuals at equilibrium $I^{*}$ (in black, bottom line) and
$\widehat{I}^{*}$ (in blue, upper line) versus ${\widehat{R}}_{0}$ and
$R_{0}$, that both appear in the horizontal axis.

Figure 4: Left panel: bifurcation diagram of the number of infected
individuals at equilibrium of the classical SIS model (2, in blue, top curve)
with transmission rate $\beta=b\beta_{S}/(a+b)+a\beta_{W}/(a+b)$ and recovery
rate $\gamma=b\gamma_{S}/(a+b)+a\gamma_{W}/(a+b)$ and system (1) in black,
bottom curve. The bifurcation parameter is $\hat{R}_{0}$ and $R_{0}$,
respectively, for
$N=12,b=35.4,\beta_{S}=1,\beta_{W}=9.1,\gamma_{S}=8.4,\gamma_{W}=15.3$ and
$a\in[0.1,10]$. Note that expression ((17)) and the parameters values explain
the aparent gap in $\widehat{I}^{*}$. Central panel, the ratio
$\widehat{I}^{*}/I^{*}$. Right panel: the ratio $\widehat{R}_{0}/R_{0}$.

###  4.3 Sensitivity analysis, sensitivity indices and epidemic control

Finally, we accomplish a sensitivity analysis of the outcome of the model to
the parameters of the model.

In Section 3 we have shown that the long term behavior of the model can fully
be described in terms of $R_{0}$ and $\Delta$. We first examine the
expressions of $R_{0}$ and $\Delta$. Next, we calculate the corresponding
sensitivity indices [6]. Then we have drawn concussion useful for control
purposes.

Grouping terms in the expression of $\Delta$ (see equation (13)) yields

$\begin{array}[]{rl}\Delta&=\left(\dfrac{\gamma_{S}}{\beta_{S}}+\dfrac{\gamma_{W}}{\beta_{W}}\right)\dfrac{1}{N}+\left(\dfrac{a}{\beta_{S}}+\dfrac{b}{\beta_{W}}\right)\dfrac{1}{N}\\\
&\\\ &=\underbrace{\dfrac{1}{R_{0,S}}+\dfrac{1}{R_{0,W}}}_{\text{Block
1}}+\underbrace{\left(\dfrac{a}{\beta_{S}}+\dfrac{b}{\beta_{W}}\right)\dfrac{1}{N}}_{\text{Block
2}}\end{array}$ (18)

Interestingly, the first block includes the basic reproduction numbers
corresponding to either no weakened class (all individuals are susceptible) or
no susceptible class (all individuals are weakened, i.e., the analogous case
with other values for $\gamma$ and $\beta$). In contrast, the second block
includes the ratio of the rates at which individuals leave the susceptible
class ($a/\beta_{S}$) or the weakened class ($b/\beta_{W}$). Two conclusions
can be drawn from (18): on the one hand, $\Delta$ depends linearly on $a$ and
$b$ (see the top right panel of Figure 5). Therefore, $\Delta$ changes
linearly with these parameters. On the other hand, $\Delta$ depends non-
linearly on the corresponding basic reproduction numbers or, ultimately, on
the transmission rates (see the bottom right panel of Figure 5).

Figure 5: Left (right, respt.) column, $R_{0}$ ($\Delta$, respt.) as function
of coefficients $a,\,b$ (top) and $\beta_{W},\,\beta_{S}$ (bottom). All the
panels display also the threshold plane $R_{0}=1$ (left column) and $\Delta=1$
(right column).

Concerning $R_{0}$, it depends linearly on $\beta_{S}$ and $\beta_{W}$ (see
expression (8) and the bottom left panel of Figure 5). On the contrary,
$R_{0}$ depends non-linearly on $a$ and $b$, although the effect is quasi-
linear (see expression (8) and the top left panel of Figure 5).

####  4.3.1 Sensitivity indices

In avoiding the communicable disease becoming endemic (or promoting
endemicity) we must control the values of $R_{0}$ and $\Delta$ (see Figure 3).
It is useful to know the relative importance of the parameters involved in the
expressions of $R_{0}$ and $\Delta$, so we can choose which of them must be
changed when developing intervention strategies.

The normalized forward sensitivity index of a variable, $u$, that depends
differentiably on a parameter $p$ is defined in [6] as:

$\gamma_{p}^{u}:=\frac{\partial u}{\partial
p}\times\frac{p}{u}=\frac{\frac{\partial u}{\partial p}}{\frac{u}{p}}$ (19)

For the convenience of the reader we sketch the derivation of the sensitivity
index (19) in Appendix 5. Its interpretation is as follows. When increasing
(or decreasing) $p_{0}$ by $\epsilon$ (meaning increasing $p_{0}$ to
$p_{0}+\epsilon p_{0}=(1+\epsilon)p_{0}$) $u$ increases (or decreases) by
$\epsilon\left.\gamma_{p}^{u}\right|_{p=p_{0}}$ times $u(p_{0})$. Let us
underline that this interpretation is local and approximated in the same sense
the Taylor’s expansion is so (see Appendix (5)).

We have used (19) to derive the analytical expression for the sensitivity
index of $R_{0}$ and $\Delta$, defined by (8) and (13) respectively, to each
of the six parameters considered in our model. We show the corresponding
expressions in Table 1.

| $R_{0}$ | $\Delta$
---|---|---
$a$ | $\frac{ab}{a+b}\frac{\beta_{w}-\beta_{s}}{b\beta_{s}+a\beta_{w}}$ | $\frac{a\beta_{w}}{(a+\gamma_{s})\beta_{w}+(b+\gamma_{w})\beta_{s}}$
$b$ | $\frac{ab}{a+b}\frac{\beta_{s}-\beta_{w}}{b\beta_{s}+a\beta_{w}}$ | $\frac{b\beta_{s}}{(a+\gamma_{s})\beta_{w}+(b+\gamma_{w})\beta_{s}}$
$\beta_{s}$ | $\frac{b\beta_{s}}{b\beta_{s}+a\beta_{w}}$ | $-\frac{(a+\gamma_{s})\beta_{w}}{(a+\gamma_{s})\beta_{w}+(b+\gamma_{w})\beta_{s}}$
$\beta_{w}$ | $\frac{a\beta_{w}}{b\beta_{s}+a\beta_{w}}$ | $-\frac{(b+\gamma_{w})\beta_{s}}{(a+\gamma_{s})\beta_{w}+(b+\gamma_{w})\beta_{s}}$
$\gamma_{s}$ | $-\frac{\gamma_{s}}{\gamma_{s}+\gamma_{w}}$ | $\frac{\gamma_{s}\beta_{w}}{(a+\gamma_{s})\beta_{w}+(b+\gamma_{w})\beta_{s}}$
$\gamma_{w}$ | $-\frac{\gamma_{w}}{\gamma_{s}+\gamma_{w}}$ | $\frac{\gamma_{w}\beta_{s}}{(a+\gamma_{s})\beta_{w}+(b+\gamma_{w})\beta_{s}}$
Table 1: Sensitivity indices of $R_{0}$ (8) and $\Delta$ (13) to recovery and
transmission rates of the NCD and the transmissible disease considered in our
model.

Assume now that $u$ depends on parameters $p_{1},\cdots,p_{n}$. Without lost
of generality, we assume that all the parameters vary simultaneously by
$\epsilon$. A direct application of the generalized Taylor’s expansion yields
that the corresponding sensitivity index is

$\sum_{i=1}^{n}\gamma_{p_{i}}^{u}$ (20)

That is to say that the sensitivity indices are additive, but (or and) the
sign of each sensitivity index matters. Note that when the sum of the
different indices is 1 (or -1) the variation caused in the variable $u$ is the
same variation introduced in the parameters times the corresponding
sensitivity index.

Most of the expressions of these sensitivity indices are complex so it is not
possible to set an order from most sensitive to least sensitive without
evaluating them at some baseline parameter values. However, some general
conclusions can be drawn.

####  4.3.2 Implications for managing disease outbreaks

Epidemiologists look at $R_{0}$ at the beginning of epidemic outbreaks [3],
[16] (but see also [8] for a less theoretical approach). $R_{0}$ depends on
the parameters of the model and a key question is that of ascertaining which
coefficients modify to get the larger change in $R_{0}$ with minimum effort.

We next derive information from the expression of the sensitivity indices
useful for disease managers (see Table 1) that may be useful for disease
managers. All the mathematical relations follow straightforward from the
expressions gathered in Table 1. We assume that
$a,\,b,\,\beta_{s},\,\beta_{w},\,\gamma_{s},\,\gamma_{w}>0$.

Bounds for the sensitivity indices On the one hand, most of the expressions
for the sensitivity indices are fractions in which the numerator is one of the
summands of the denominator. As all the parameters are positive quantities,
the absolute value of these indices is less than one. More specifically, Table
2 shows infimum and supremum values of sensitivity indices of $R_{0}$ and
$\Delta$.

| Sensitivity indices of $R_{0}$ | Sensitivity indices of $\Delta$
---|---|---
| Infimum value | Supremum value | Infimum value | Supremum Value
$a$ | - | - | 0 | 1
$b$ | - | - | 0 | 1
$\beta_{s}$ | 0 | 1 | -1 | 0
$\beta_{w}$ | 0 | 1 | -1 | 0
$\gamma_{s}$ | -1 | 0 | 0 | 1
$\gamma_{s}$ | -1 | 0 | 0 | 1
Table 2: Infimum and supremum values of sensitivity index of $R_{0}$ and
$\Delta$. The cells are empty when no meaningful bounds can be provided.

As shown in Table 2, the supremum value of the sensitivity indices is $1$. So,
it is not possible to change $R_{0}$ (or $\Delta$) by an amount bigger than
the change $\epsilon$ introduced in the parameter times $R_{0}$ (or
$\Delta$).The interpretation of the infimum value $-1$ is the same, but in
that case, the modification introduced in the parameters and the change in
$R_{0}$ (or $\Delta$) have the opposite direction. In addition, note that we
are considering a Taylor expansion to first order, so no big changes can be
targeted.

Comparing communicable disease management strategies Tables 1 and 2 allow us
to compute how much modifying a single coefficient of the model makes $R_{0}$
or $\Delta$ vary. However, it is possible to act on more than one coefficient
of the model at once. Table 3 gathers combinations of parameters to be
modified simultaneously and equally that produce minimal, none, or maximal
responses on $R_{0}$ and $\Delta$.

$R_{0}$ | $\Delta$
---|---
$\gamma_{\beta_{s}}^{R_{0}}+\gamma_{\beta_{w}}^{R_{0}}=1$ | $\gamma_{\beta_{s}}^{\Delta}+\gamma_{\beta_{w}}^{\Delta}=-1$
$\gamma_{\gamma_{s}}^{R_{0}}+\gamma_{\gamma_{w}}^{R_{0}}=-1$ | –
$\gamma_{\gamma_{s}}^{\Delta}+\gamma_{\gamma_{w}}^{\Delta}+\gamma_{a}^{\Delta}+\gamma_{b}^{\Delta}=-1$ | $\gamma_{\gamma_{s}}^{\Delta}+\gamma_{\gamma_{w}}^{\Delta}+\gamma_{a}^{\Delta}+\gamma_{b}^{\Delta}=1$
$\gamma_{a}^{R_{0}}+\gamma_{b}^{R_{0}}=0$ | –
$\gamma_{\beta_{s}}^{R_{0}}+\gamma_{\beta_{w}}^{R_{0}}+\gamma_{\gamma_{s}}^{R_{0}}+\gamma_{\gamma_{w}}^{R_{0}}+\gamma_{a}^{R_{0}}+\gamma_{b}^{R_{0}}=0$ | $\gamma_{\beta_{s}}^{\Delta}+\gamma_{\beta_{w}}^{\Delta}+\gamma_{\gamma_{s}}^{\Delta}+\gamma_{\gamma_{w}}^{\Delta}+\gamma_{a}^{\Delta}+\gamma_{b}^{\Delta}=0$
Table 3: Expressions of sums of indices when different combinations of
parameters are simultaneously modified producing minimal, none or maximal
response on $R_{0}$ and $\Delta$.

We next examine how the expressions gathered in Table 3 can be used to decide
on epidemics management strategies.

A key question from the point of view of managing a CD that of deciding to act
either on a target population (for instance, weakened individuals) or equally
on all the susceptible individuals regardless of their status. More
specifically, the question is: What will produce a larger change in $R_{0}$,
an effect $\epsilon_{1}$ applied only on (say) $\beta_{w}$, or a weaker effect
$\epsilon_{2}<\epsilon_{1}$ applied on both $\beta_{s}$ and $\beta_{w}$? This
question is equivalent to compare
$\gamma^{R_{0}}_{\beta_{w}}\epsilon_{1}R_{0}$ to
$(\gamma^{R_{0}}_{\beta_{w}}+\gamma^{R_{0}}_{\beta_{s}})\epsilon_{2}R_{0}$.
Note that

$\gamma_{\beta_{s}}^{R_{0}}+\gamma_{\beta_{w}}^{R_{0}}=1.$ (21)

Direct calculations yield that

$\gamma^{R_{0}}_{\beta_{w}}\epsilon_{1}R_{0}<(\gamma^{R_{0}}_{\beta_{w}}+\gamma^{R_{0}}_{\beta_{s}})\epsilon_{2}R_{0}\qquad\Leftrightarrow\qquad\frac{\epsilon_{1}}{\epsilon_{2}}<\frac{1}{\gamma^{R_{0}}_{\beta_{w}}}$
(22)

provided (21). Note that $\epsilon$ must be negative in order to reduce
transmission.

Analogous questions can be addressed related to those expressions summing up
to $-1$ or $0$.

Are unexpected management effects possible? We already know that the endemic
states bifurcate from the disease-free scenario as $R_{0}$ crosses the
threshold value $1$. The bifurcation can be either subcritical ($R_{0}<1$ and
$\Delta<1$) or supercritical ($R_{0}\geq 1$ and $\Delta\geq 1$). The
bifurcation direction would make a huge difference, as in the subcritical case
$R_{0}<1$ does not lead necessarily to a disease-free scenario. Changing any
parameter of the model will make vary simultaneously $R_{0}$ and $\Delta$
which would result in an unexpected outcome.

For instance, let us assume that the conditions are such that $R_{0}>1$ and
$\Delta>1$, that is, the system is in the endemic disease scenario. Assume
also that efforts are put into modifying the value of some coefficients to
push $R_{0}$ below $1$. Then, as a result, can $\Delta$ also go below $1$,
undergoing a subcritical bifurcation? It would happen if reducing $R_{0}$ (by
any means) would entail a simultaneous reduction in $\Delta$. A necessary
condition is that the corresponding sensitivity indices have the same sing. It
is apparent that it is not possible if the efforts are put in modifying
coefficients $\beta_{s}$, $\beta_{w}$, $\gamma_{s}$ or $\gamma_{w}$ (see Table
1).

However, it is also apparent that the sign of the sensitivity indices of
$R_{0}$ and $\Delta$ with respect to $a$ (respectively $b$) is the same
provided $\beta_{w}>\beta_{s}$ (respectively, if $\beta_{w}<\beta_{s}$).
Indeed direct calculations show that $R_{0}>1$ and $\Delta>1$ for $N=70$,
$a=0.27$, $b=15$, $\beta_{S}=3.02$, $\beta_{W}=7.2$, $\gamma_{S}=20$,
$\gamma_{W}=6.9$ but $R_{0}<1$ and $\Delta<1$ keeping all the previous
parameter values but $a=0.17$.

We hope the results presented herein would promote further research. On the
one hand, we hope experimental scientists find this research interesting and
would test the model at their laboratory. On the other hand, more realistic
extensions of system (1) should be of interest.

## Declarations:

Funding: The work of E. Venturino has been partially supported by the projects
“Metodi numerici in teoria delle popolazioni”, “Metodi numerici nelle scienze
applicate” of the Dipartimento di Matematica “Giuseppe Peano” of the
Università di Torino and the program “Giner de los Ríos” of the Universidad de
Alcalá. The work of M. Marvá has been partially supported by Ministerio de
Economía y Competitividad (Spain), project MTM2014-56022-C2-1 and Ministerio
de Ciencia, Innovación y Universidades (Spain), project
RTI2018-096884-B-C32-P.

Conflicts of interest/Competing interests: the authors of this work certify
that they have NO affiliations with or involvement in any organization or
entity with any financial interest or non-financial interest in the subject
matter or materials discussed in this manuscript.

## References

  * [1] P. Auger, R. Bravo de la Parra, J.C. Poggiale, E. Sánchez, and L. Sanz.  Aggregation methods in dynamical systems and applications in population and community dynamics.  Physics of Life Reviews, 5(2):79–105, jun 2008.
  * [2] M. Begon, M. Bennett, R. G. Bowers, N. P. French, S. M. Hazel, and J. Turner.  A clarification of transmission terms in host-microparasite models: numbers, densities and areas.  Epidemiology and Infection, 129(1):147–153, aug 2002.
  * [3] S. Boonpatcharanon, J. Heffernan, and H. Jankowski.  Estimating the basic reproduction number at the beginning of an outbreak under incomplete data.  Plos One, jul 2021.
  * [4] F. Brauer.  Backward bifurcations in simple vaccination models.  Journal of Mathematical Analysis and Applications, 298(2):418–431, oct 2004.
  * [5] B. Buonomo.  A note on the direction of the transcritical bifurcation in epidemic models.  Nonlinear Analysis: Modelling and Control, 20(1):38–55, mar 2015\.
  * [6] N. Chitnis, J.M. Hyman, and J.M. Cushing.  Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model.  Bulletin of Mathematical Biology, 70(5):1272–1296, feb 2008.
  * [7] GBD 2015 Risk Factors Collaborators.  Global, regional, and national comparative risk assessment of 79 behavioural, environmental and occupational, and metabolic risks or clusters of risks, 1990–2015: a systematic analysis for the global burden of disease study 2015.  Lancet, 388(10053):1659–1724, 2016.
  * [8] Paul L. Delamater, Erica J. Street, Timothy F. Leslie, Y. Tony Yang, and Kathryn H. Jacobsen.  Complexity of the basic reproduction number ($r_{0}$).  Emerging Infectious Diseases, 25(1):1–4, jan 2019.
  * [9] M.M. Coates et al. Burden of non-communicable diseases from infectious causes in 2017: a modelling study.  The Lancet Global Health, 8(12):e1489–e1498, dec 2020.
  * [10] Z. Feng, C. Castillo-Chavez., and A.F. Capurro.  A model for tuberculosis with exogenous reinfection.  Theoretical Population Biology, 57(3):235–247, may 2000.
  * [11] E. Garde, B. Ibrahim, A.T. Kovács, and S. Schuster.  Differential equation-based minimal model describing metabolic oscillations in bacillus subtilis biofilms.  Royal Society Open Science, 7(2):190810, feb 2020.
  * [12] D. Greenhalgh, O. Diekmann, and M.C.M. de Jong.  Subcritical endemic steady states in mathematical models for animal infections with incomplete immunity.  Mathematical Biosciences, 165(1):1–25, may 2000.
  * [13] A.B. Gumel.  Causes of backward bifurcations in some epidemiological models.  Journal of Mathematical Analysis and Applications, 395(1):355–365, nov 2012.
  * [14] A.B. Gumel, E.A. Iboi, C.N. Ngonghala, and E.H. Elbasha.  A primer on using mathematics to understand COVID-19 dynamics: Modeling, analysis and simulations.  Infectious Disease Modelling, 6:148–168, 2021.
  * [15] K.P. Hadeler and P. Van den Driessche.  Backward bifurcation in epidemic control.  Mathematical Biosciences, 146(1):15–35, nov 1997.
  * [16] M.J. Keeling and P. Rohani.  Modeling Infectious Diseases in Humans and Animals.  PRINCETON UNIV PR, October 2007.
  * [17] K. Mark and T.D. Dobromir.  The dynamics of a simple, risk-structured HIV model.  Mathematical Biosciences and Engineering, 17(4):4184–4209, 2020\.
  * [18] M. Martcheva.  An Introduction to Mathematical Epidemiology.  Springer-Verlag GmbH, October 2015.
  * [19] M. Marvá, R. Bravo de la Parra, and E. Venturino.  Modelling the role of opportunistic diseases in coinfection.  Mathematical Modelling of Natural Phenomena, 13(3):28, 2018.
  * [20] M. Marvá, E. Venturino, and R. Bravo de la Parra.  A time scales approach to coinfection by opportunistic diseases.  Journal of Applied Mathematics, 2015:1–10, 2015.
  * [21] N.K.-D. Ohene Opoku, G. Bader, and E. Fiatsonu.  Controlling crime with its associated cost during festive periods using mathematical techniques.  Chaos, Solitons $\&$ Fractals, 145:110801, apr 2021.
  * [22] A. Prüss-Ustün, E. van Deventer, P. P. Mudu, D. Campbell-Lendrum, C. Vickers, I Ivanov, F. Forastiere, Gumy S., Dora C., H Adair-Rohani, and M. Neira.  Environmental risks and non-communicable diseases.  BMJ, page l265, jan 2019.
  * [23] S. M. Raimundo, H.M. Yang, and E. Massad.  Contagious criminal career models showing backward bifurcations: Implications for crime control policies.  Journal of Applied Mathematics, 2018:1–16, aug 2018.
  * [24] T.C. Reluga, J. Medlock, and A.S. Perelson.  Backward bifurcations and multiple equilibria in epidemic models with structured immunity.  Journal of Theoretical Biology, 252(1):155–165, may 2008.
  * [25] E. Venturino and S. Pandey.  A TB model: Is disease eradication possible in india?  Mathematical Biosciences and Engineering, 15(1):233–254, may 2017\.
  * [26] WHO.  The top 10 causes of death.  World Health Organization., Retrieved June 2021.
  * [27] E.B. Wong and et al.  Convergence of infectious and non-communicable disease epidemics in rural south africa: a cross-sectional, population-based multimorbidity study.  The Lancet Global Health, 9(7):e967–e976, jul 2021.

##  5 Appendix

A first approach to the idea of the sensitivity of $u$ to $p$ is using the
derivative of $u$ with respect to $p$. However, doing so does not allow to
fairly compare the sensitivity of $p$ to two different parameters if those
parameters are expressed in different units. Defining the sensitivity index as
in (19) fixes this problem (see the most right hand side expression).
Furthermore, consider that $u$ depends on the parameters $p_{1},\cdots,p_{n}$.
We may assume without lost of generality that $p_{2},\cdots,p_{n}$ are held
constant, that is equivalent to assume that $u$ depends only on $p=p_{1}$. The
Taylor’s expansion approximation of $u(p)$ to the first order at $p=p_{0},$ is
given by:

$u(p)\approx u(p_{0})+\left.\frac{\partial u(p)}{\partial
p}\right|_{p_{0}}(p-p_{0})$ (23)

When varying the parameter $p$ by an amount of $\epsilon=(p-p_{0})/p_{0}$,
that is to say, from $p_{0}$ to $p_{0}+\epsilon p_{0}$ (23) becomes

$u(p_{0}+\epsilon p_{0})-u(p_{0})\approx\left.\frac{\partial u(p)}{\partial
p}\right|_{p_{0}}\epsilon p_{0}.$ (24)

Taking into account the definition of the sensitivity index (19), multiplying
and dividing the right hand side of (24) by $u(p_{0})$ yields:

$u(p_{0}+\epsilon
p_{0})-u(p_{0})\approx\left.\gamma_{p}^{u}\right|_{p_{0}}\epsilon u(p_{0})$
(25)