PLoS ONE
Determination of critical community size from an HIV/AIDS model
DOI 10.1371/journal.pone.0244543 , Volume: 16 , Issue: 1
Article Type: research-article, Article History
•
•
• Altmetric

### Notes

Abstract

After an epidemic outbreak, the infection persists in a community long enough to engulf the entire susceptible population. Local extinction of the disease could be possible if the susceptible population gets depleted. In large communities, the tendency of eventual damp down of recurrent epidemics is balanced by random variability. But, in small communities, the infection would die out when the number of susceptible falls below a certain threshold. Critical community size (CCS) is considered to be the mentioned threshold, at which the infection is as likely as not to die out after a major epidemic for small communities unless reintroduced from outside. The determination of CCS could aid in devising systematic control strategies to eradicate the infectious disease from small communities. In this article, we have come up with a simplified computation based approach to deduce the CCS of HIV disease dynamics. We consider a deterministic HIV model proposed by Silva and Torres, and following Nåsell, introduce stochasticity in the model through time-varying population sizes of different compartments. Besides, Metcalf’s group observed that the relative risk of extinction of some infections on islands is almost double that in the mainlands i.e. infections cease to exist at a significantly higher rate in islands compared to the mainlands. They attributed this phenomenon to the greater recolonization in the mainlands. Interestingly, the application of our method on demographic facts and figures of countries in the AIDS belt of Africa led us to expect that existing control measures and isolated locations would assist in temporary eradication of HIV infection much faster. For example, our method suggests that through systematic control strategies, after 7.36 years HIV epidemics will temporarily be eradicated from different communes of island nation Madagascar, where the population size falls below its CCS value, unless the disease is reintroduced from outside.

Das, Ghosh, Banerjee, Pyne, Chattopadhyay, Mukhopadhyay, and Forgoston: Determination of critical community size from an HIV/AIDS model

## Introduction

Periodicity in the recurrence of local epidemics has attracted mathematicians and epidemiologists for centuries. Soper [1] pointed out that the time to extinction of infection and the size of the susceptible population are possibly the responsible candidates behind the periodicity. The persistence of a disease in a community, often measured by the time for an epidemic to die out, is a major concern to date. After an epidemic outbreak, the infection persists in the community long enough to engulf the entire susceptible population. Local extinction of the disease could be possible if the susceptible population gets depleted. The success of Hamer-Soper model [1] in predicting periodicity of epidemics from simple assumptions, led Bartlett [24] to observe a threshold for the susceptible population size. If the susceptible population is equal to that threshold, there is a 50% chance that the infection would cease to exist after the corresponding mean fade-out time unless the disease is reintroduced from outside. He analyzed this threshold using a semi-stochastic version of the Hamer-Soper model and established time to extinction of the epidemic using quasi stationarity. Bartlett’s model sufficiently explained the undamped epidemic outbreaks and its extinction. Bartlett applied the semi-stochastic model for the recurrent measles outbreak in England and Wales. He observed that in large communities, the tendency of eventual damp down of recurrent epidemics is balanced by random variability; but, in small communities, the infection would die out when the size of the susceptible population falls below a certain threshold. These led Bartlett to define Critical Community Size (CCS) as a threshold size of the susceptible population in a small community, at which infection is as likely as not to die out, after a major epidemic, unless reintroduced from outside. The idea is that, if the size of susceptible population drops below that threshold, the pathogen might cease to exist in the population.

Apart from the works by Andersson and Britton [5], scarcely significant mathematical contributions were made in this direction until Nåsell [6] reformulated CCS for Bartlett’s Susceptible-Infected (SI) model using quasi stationarity under more realistic assumptions. Nåsell [6, 7] used a fully stochastic model to study the time to extinction of an epidemic for a number of disease models. The model and method proposed by Nåsell [6] is based on a standard SI model that involves a system of two differential equations. This model can easily be extended to an SIR model where the population size remains constant at any time, as this SIR model can be reduced to a two-equation model and the above method is applicable without any significant modification. Nåsell [6] also focused to unleash the potential of CCS as a control strategy for the epidemic outbreak. Later, Metcalf et al. [8] mentioned this potential of CCS. But, calculation of CCS is limited mostly to models based on a two-equation system of differential equations or measles type diseases. One predominant reason is that the analytical formulation of CCS for complex infectious diseases is quite cumbersome. The absence of any method to calculate CCS for higher-order models motivated us to develop a cost (in terms of time) effective technique to evaluate CCS. In this paper, we propose a simplified computation based approach for rapid calculation of CCS, that could emphasize prompt formulation of control strategies and/or policies at the time of disease outbreak.

Technically, this paper determines the time to extinction of an epidemic outbreak theoretically and also develops a computation based approximate formulation of CCS using the diffusion approximation. The application of our method is free from any restriction on the number of variables in the system that defines the disease dynamics. To demonstrate our method we consider an HIV model by Silva, Torres, and Djordjevic [911] that describes the disease dynamics using a system of four differential equations. Since infectious diseases mostly have a stochastic nature [1214], we consider a stochastic version of the HIV disease dynamics. Following Nåsell, we introduce stochasticity through time-varying population sizes of different compartments in the deterministic HIV model [9]. Next, we calculate the mean time to extinction of the HIV epidemic and compute approximate value of CCS using diffusion approximation based on the Ornstein-Uhlenbeck process.

In this paper, we have cited the HIV epidemic case of Madagascar, in Sub-Saharan Africa because of a special reason. In spite of less remarkable progress in the 90-90-90 HIV treatment target of the Joint United Nations Programme on HIV and AIDS (UNAIDS), Madagascar has shown a much lower prevalence and incidence rate of HIV infection, compared to some other African countries. It is interesting to note that, despite being in the part of Africa where HIV is in a state of epidemic, it has much lower rates of prevalence and incidence. Incidentally, Madagascar is situated in a geographically isolated location, separated by 60 km of water from the other eastern and southern countries that lie in the mainland of the African continent. This observation echoes a fact from Metcalf et al. [8], that says the relative risk of extinction of some infectious diseases on islands is almost double that of those in the mainland. They also suggest that persistence of the infection within a region in the mainland is fuelled constantly by recolonization and infections cease to exist at significantly higher rates in islands where access to and from the mainland is much less. This fascinating phenomenon led us to postulate that, in conjunction with UNAIDS controlling strategies if it is possible to keep the HIV susceptible population of Madagascar small or rather below CCS, effective and rapid eradication of HIV might be possible there, after the mean extinction time of the disease. This is primarily based on the fact that the isolated geographical location of Madagascar might lead to lesser chances of reintroduction of the infection.

We, therefore, applied our method on the demographic facts and figures of Madagascar to obtain the CCS value and mean time to extinction of the disease. Our calculations suggest that, unless the disease is reintroduced from outside, after 7.36 years HIV epidemics might fade out from different cities or communes of the island nation where the susceptible population size drops below 4585, the CCS value of Madagascar. This time period of 7.36 years could appear astonishing as HIV has become a long-standing epidemic situation in the countries of the Sub-Saharan AIDS belt. But there is a strong possibility that the idea of CCS would propel the controlling strategies of UNAIDS to eradicate the disease in the near future. This work suggests that the systematic spreading of control measures could lead to accelerated progress to achieve the 90-90-90 target of UNAIDS for HIV eradication. Such progress shall prevent the HIV epidemic in Madagascar unless the infection is reintroduced from outside. Moreover, countries like Uganda, Malawi (located in Eastern and Southern Africa) are separated from their neighboring countries by water bodies such as fresh water lakes. These countries have lower rates of prevalence and incidence compared to landlocked countries such as Eswatini (in Southern Africa). Notably, these rates are still much higher than that of Madagascar.

To date, CCS has largely been a part of mathematical theory only. But implementing CCS as controlling strategies of epidemics might prove beneficial in the future. Therefore, to utilize the potential that CCS holds, it is important to determine its value promptly for any disease dynamics. This motivated us to calculate CCS (approximately) using simplified computations. Our method will help in finding country-specific CCS and time to extinction of the disease, using real demographic facts and figures for HIV infection in different countries as parameter values. We emphasize the potential utility of CCS to guide public health practitioners to develop appropriate area-specific control strategies and interventions especially when they are limited, before or after an epidemic outbreak.

## Materials and methods

### Model formulation

We work with the model proposed by Silva and Torres [9] that explains the dynamics of HIV infection assuming homogenous mixing within communities. Here the entire population is divided into four compartments. These compartments are mutually exclusive in the sense that no person can belong to more than one compartment at any time point. The four compartments are: susceptible individuals (S), individuals infected with HIV but having no clinical symptoms of AIDS (I), individuals infected with AIDS and having ART treatment as they are supposed to belong to chronic stage (C), and individuals infected with HIV having clinical symptoms of AIDS (A). Note that a person belonging to I may transmit HIV to other individuals although he/she has no symptoms of AIDS. Under this situation, the model we consider is,

$\frac{dS}{dt}=\Lambda -\beta \left(I\left(t\right)+{\eta }_{C}C\left(t\right)+{\eta }_{A}A\left(t\right)\right)\frac{S\left(t\right)}{N}-\mu S\left(t\right)$
$\frac{dI}{dt}=\beta \left(I\left(t\right)+{\eta }_{C}C\left(t\right)+{\eta }_{A}A\left(t\right)\right)\frac{S\left(t\right)}{N}-\left(\rho +\varphi +\mu \right)I\left(t\right)+\omega C\left(t\right)+\alpha A\left(t\right)$
$\frac{dC}{dt}=\varphi I\left(t\right)-\left(\omega +\mu \right)C\left(t\right)$
$\frac{dA}{dt}=\rho I\left(t\right)-\left(\alpha +\mu +d\right)A\left(t\right)$

Here β represents the contact rate for HIV transmission, ηA (≥ 1) is the relative infectiousness of individuals with AIDS symptoms compared to those infected with HIV but no AIDS symptoms, and ηC (≤ 1) is the rate of partial restoration of immune function of HIV infected individuals who are correctly treated under ART. ρ and ϕ are the rates at which individuals in class I move to A and C respectively. α and ω are the respective rates at which individuals from class A and class C move to I. We assume that d is the HIV induced death rate while μ is the natural death rate in the population.

We consider the total population as N(t) = S(t) + I(t) + A(t) + C(t) implying that $\frac{dN\left(t\right)}{dt}=\Lambda -\mu N\left(t\right)-dA\left(t\right)$. It can be easily shown that $N\left(t\right)\le \frac{\Lambda }{\mu }$.

### Basic reproduction number

We use the approach of next generation matrix [15] to calculate R0. Note that the rate of new infections coming in the compartments I, C and A are $\beta \left(I+{\eta }_{C}C+{\eta }_{A}A\right)\frac{S}{N}$, 0 and 0 respectively. On the other hand, new infections moving out of the compartments I, C, and A are (ρ + ϕ + μ)IωCαA, (ω + μ)CϕI and (α + μ + d)AρI respectively. Thus following [9, 15], we have basic reproduction number (R0) as:

$\begin{array}{c}\hfill {R}_{0}=\frac{{S}^{0}\beta }{N}\frac{{\xi }_{2}\left({\xi }_{1}+\rho {\eta }_{A}\right)+{\eta }_{C}\varphi {\xi }_{1}}{\mu \left({\xi }_{2}+\left(\rho +{\xi }_{1}\right)+\varphi {\xi }_{1}+\rho d\right)+\rho \omega d}=\frac{{S}^{0}}{N}\frac{\zeta }{\gamma }\approx \frac{\zeta }{\gamma }\end{array}$
where ζ = β(ξ2(ξ1 + ρηA) + ηC ϕξ1), γ = μ(ξ2 + (ρ + ξ1) + ϕξ1 + ρd) + ρωd, ξ1 = α + μ + d and ξ2 = ω + μ.

### Stochastic model and quasi-stationarity

Evaluation of CCS depends on the time to extinction of disease based on the given model. In order to find the time to extinction, we now construct the fully stochastic version of the model (1)—(4). Next, we write the Kolmogorov forward equation for this stochastic model. This would be the key step for further downstream analysis.

First, we note the nature of transitions and the respective transition rates from one compartment to another (Table 1).

Table 1
Chart for transition rates.
EventTransitionTransition rate
Immigration of susceptibles (S)(s, i, c, a) → (s + 1, i, c, a)λ1(s, i, c, a) = μN
Death of susceptibles (S)(s, i, c, a) → (s − 1, i, c, a)λ2(s, i, c, a) = μs
Susceptible (S) to Infected (I)(s, i, c, a) → (s − 1, i + 1, c, a)λ3(s, i, c, a) = βsi/N
Susceptible (S) to Infected under treatment (C)(s, i, c, a) → (s − 1, i, c + 1, a)λ4(s, i, c, a) = βηC sc/N
Susceptible (S) to Infected with AIDS (A)(s, i, c, a) → (s − 1, i, c, a + 1)λ5(s, i, c, a) = βηA sa/N
Infected (I) to Infected under treatment(C)(s, i, c, a) → (s, i − 1, c + 1, a)λ6(s, i, c, a) = ϕi
Infected (I) to Infected with AIDS (A)(s, i, c, a) → (s, i − 1, c, a + 1)λ7(s, i, c, a) = ρi
Death of Infected (I)(s, i, c, a) → (s, i − 1, c, a)λ8(s, i, c, a) = μi
Infected under treatment (C) to infected (I)(s, i, c, a) → (s, i + 1, c − 1, a)λ9(s, i, c, a) = ωc
Death of Infected under treatment(C)(s, i, c, a) → (s, i, c − 1, a)λ10(s, i, c, a) = μc
Infected with AIDS symptoms (A) to Infected (I)(s, i, c, a) → (s, i + 1, c, a − 1)λ11(s, i, c, a) = αa
Death of Infected with AIDS symptoms (A)(s, i, c, a) → (s, i, c, a − 1)λ12(s, i, c, a) = (μ + d)a

Here, s = S/N, i = I/N, c = C/N, a = A/N. Then the Kolmogorov forward equation becomes,

$\begin{array}{cc}{p}_{s,i,c,a}^{\prime }\left(t\right)\hfill & ={\lambda }_{1}\left(s-1,i,c,a\right){p}_{s-1,i,c,a}+{\lambda }_{2}\left(s+1,i,c,a\right){p}_{s+1,i,c,a}\hfill \\ & +{\lambda }_{3}\left(s+1,i-1,c,a\right){p}_{s+1,i-1,c,a}+{\lambda }_{4}\left(s+1,i,c-1,a\right){p}_{s+1,i,c-1,a}\hfill \\ & +{\lambda }_{5}\left(s+1,i,c,a-1\right){p}_{s+1,i,c,a-1}+{\lambda }_{6}\left(s,i+1,c-1,a\right){p}_{s,i+1,c-1,a}\hfill \\ & +{\lambda }_{7}\left(s,i+1,c,a-1\right){p}_{s,i+1,c,a-1}+{\lambda }_{8}\left(s,i+1,c,a\right){p}_{s,i+1,c,a}\hfill \\ & +{\lambda }_{9}\left(s,i-1,c+1,a\right){p}_{s,i-1,c+1,a}+{\lambda }_{10}\left(s,i,c+1,a\right){p}_{s,i,c+1,a}\hfill \\ & +{\lambda }_{11}\left(s,i-1,c,a+1\right){p}_{s,i-1,c,a+1}+{\lambda }_{12}\left(s,i,c,a+1\right){p}_{s,i,c,a+1}-\kappa \left(s,i,c,a\right){p}_{s,i,c,a}\hfill \end{array}$
where $\kappa \left(s,i,c,a\right)=\sum _{i=1}^{12}{\lambda }_{i}\left(s,i,c,a\right)$ and ps,i,c,a(t) = P[S(t) = s, I(t) = i, C(t) = c, A(t) = a].

Now, we find quasi-stationarity by conditioning on non-extinction of the disease. So we write,

$\begin{array}{c}\hfill {q}_{s,i,c,a}\left(t\right)=P\left[S\left(t\right)=s,I\left(t\right)=i,C\left(t\right)=c,A\left(t\right)=a|\left(I\left(t\right),C\left(t\right),A\left(t\right)\right)\ne \left(0,0,0\right)\right]=\frac{{p}_{s,i,c,a}\left(t\right)}{1-{p}_{•000}\left(t\right)}\end{array}$
where ${p}_{•000}\left(t\right)=\sum _{s=0}^{\infty }P\left[S\left(t\right)=s,I\left(t\right)=0,C\left(t\right)=0,A\left(t\right)=0\right]=\sum _{s=0}^{\infty }{p}_{s,0,0,0}\left(t\right)$.

Differentiating qs,i,c,a(t) with respect to t, we have,

$\begin{array}{c}\hfill {q}_{s,i,c,a}^{\prime }\left(t\right)=\frac{{p}_{s,i,c,a}^{\prime }\left(t\right)}{1-{p}_{•000}\left(t\right)}+\frac{{p}_{s,i,c,a}\left(t\right)}{{\left(1-{p}_{•000}\left(t\right)\right)}^{2}}.{p}_{•000}^{\prime }\left(t\right).\end{array}$
Now, from Eq (6), we have, after simplification,
$\begin{array}{cc}{p}_{•000}^{\prime }\left(t\right)\hfill & =\sum _{s=0}^{\infty }\left[\mu N{p}_{s-1,0,0,0}\left(t\right)+\mu \left(s+1\right){p}_{s+1,0,0,0}\left(t\right)+\mu \left(0+1\right){p}_{s,1,0,0}\left(t\right)\hfill \\ & +\mu \left(0+1\right){p}_{s,0,1,0}\left(t\right)+\left(\mu +d\right)\left(0+1\right){p}_{s,0,0,1}\left(t\right)\right]-\left(\mu N+\mu s\right)\sum _{s=0}^{\infty }{p}_{s,0,0,0}\left(t\right)\hfill \\ & =\mu {p}_{•1\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)+\mu {p}_{•0\phantom{\rule{0.166667em}{0ex}}1\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)+\left(d+\mu \right){p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}1}\left(t\right)={p}_{•}^{\left(d,\mu \right)}\left(t\right)\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\left(say\right).\hfill \end{array}$
From Eqs (7) and (8) we have,
$\begin{array}{cc}{q}_{s,i,c,a}^{\prime }\left(t\right)\hfill & =\frac{{p}_{s,i,c,a}^{\prime }\left(t\right)}{1-{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)}+\frac{{p}_{s,i,c,a}\left(t\right)}{{\left(1-{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)\right)}^{2}}.{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}^{\prime }\left(t\right)\hfill \\ & =\frac{{p}_{s,i,c,a}^{\prime }\left(t\right)}{1-{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)}+\frac{{p}_{s,i,c,a}\left(t\right)}{1-{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)}.\frac{{p}_{•}^{\left(d,\mu \right)}\left(t\right)}{1-{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)}\hfill \\ & =\frac{{p}_{s,i,c,a}^{\prime }\left(t\right)}{1-{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)}+\frac{{p}_{s,i,c,a}\left(t\right)}{1-{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)}.{q}_{•}^{\left(d,\mu \right)}\left(t\right),\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\left(say\right)\hfill \end{array}$
$\begin{array}{c}\text{where}\phantom{\rule{4pt}{0ex}}{q}_{•}^{\left(d,\mu \right)}\left(t\right)=\frac{{p}_{•}^{\left(d,\mu \right)}\left(t\right)}{1-{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)}\hfill \end{array}$
Now, following Nåsell [6, 7] we have,
$\begin{array}{cc}& {q}_{s,i,c,a}^{\prime }\left(t\right)=0\hfill \\ & ⇒{p}_{s,i,c,a}^{\prime }\left(t\right)=-\frac{{p}_{s,i,c,a}\left(t\right)}{\left(1-{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)\right)}.{q}_{•}^{\left(d,\mu \right)}\left(t\right)\left(1-{p}_{•0\phantom{\rule{0.166667em}{0ex}}0\phantom{\rule{0.166667em}{0ex}}0}\left(t\right)\right)=-{q}_{•}^{\left(d,\mu \right)}\left(t\right){p}_{s,i,c,a}\left(t\right)\hfill \\ & \therefore {p}_{s,i,c,a}\left(t\right)=c{e}^{-{q}_{•}^{\left(d,\mu \right)}.t}={q}_{•}^{\left(d,\mu \right)}{e}^{-{q}_{•}^{\left(d,\mu \right)}.t}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\hfill \end{array}$

Let τQ be the time to extinction when the initial distribution is equal to the quasi-stationary distribution. Therefore, we have,

$\begin{array}{c}\hfill E\left({\tau }_{Q}\right)=\frac{1}{{q}_{•}^{\left(d,\mu \right)}}.\end{array}$

### Equilibrium points

The disease-free equilibrium is obtained as: ${\Sigma }_{0}=\left({S}^{0},{I}^{0},{C}^{0},{A}^{0}\right)=\left(\frac{\Lambda }{\mu },0,0,0\right)$. To find the other endemic equilibrium, if exists, we put N = N(0), x1(t) = S(t)/N, x2(t) = I(t)/N, x3(t) = C(t)/N, and x4(t) = A(t)/N . The equilibrium point is obtained by equating the first differentiation in Eqs (1)–(4) to zero [9, 10], i.e.

$\begin{array}{ccc}& {x}_{1}^{\prime }\left(t\right)=\mu -\beta \left({x}_{2}\left(t\right)+{\eta }_{C}{x}_{3}\left(t\right)+{\eta }_{A}{x}_{4}\left(t\right)\right){x}_{1}-\mu {x}_{1}\left(t\right)=0& \end{array}$
$\begin{array}{ccc}& {x}_{2}^{\prime }\left(t\right)=\beta \left({x}_{2}\left(t\right)+{\eta }_{C}{x}_{3}\left(t\right)+{\eta }_{A}{x}_{4}\left(t\right)\right){x}_{1}-\left(\rho +\varphi +\mu \right){x}_{2}\left(t\right)+\omega {x}_{3}\left(t\right)+\alpha {x}_{4}\left(t\right)=0& \end{array}$
$\begin{array}{ccc}& {x}_{3}^{\prime }\left(t\right)=\varphi {x}_{2}\left(t\right)-\left(\omega +\mu \right){x}_{3}\left(t\right)=0& \end{array}$
$\begin{array}{cc}& {x}_{4}^{\prime }\left(t\right)=\rho {x}_{2}\left(t\right)-\left(\alpha +\mu +d\right){x}_{4}\left(t\right)=0\end{array}$

For simplicity we use the notations: xj(t) = xj for j = 1, …, 4, ξ1 = α + μ + d, ξ2 = ω + μ, γ = ξ1 ξ2(ρ + ϕ + μ) − ωϕξ1ραξ2, ζ = β(ξ1 ξ2+ ηC ϕξ1+ ηA ρξ2).

Then solving Eqs (13)–(16), we have the endemic equilibrium as:

${\stackrel{^}{x}}_{1}=\frac{\gamma }{\zeta }=\frac{1}{{R}_{0}}$
${\stackrel{^}{x}}_{2}=\frac{\left(\zeta -\gamma \right)}{\zeta \gamma }.\mu {\xi }_{1}{\xi }_{2}=\left({R}_{0}-1\right)\frac{\mu }{\zeta }.{\xi }_{1}{\xi }_{2}$
${\stackrel{^}{x}}_{3}=\frac{\left(\zeta -\gamma \right)}{\zeta \gamma }.\mu \varphi {\xi }_{1}=\left({R}_{0}-1\right)\frac{\mu }{\zeta }.\varphi {\xi }_{1}$
${\stackrel{^}{x}}_{4}=\frac{\left(\zeta -\gamma \right)}{\zeta \gamma }.\mu \rho {\xi }_{2}=\left({R}_{0}-1\right)\frac{\mu }{\zeta }.\rho {\xi }_{2}$

### Diffusion approximation

We consider a diffusion approximation to the stochastic version of our model (1)—(4). This would allow us to approximately evaluate the quasi-stationary distribution using a multivariate normal distribution when N is large and R0 is greater than 1. First we find the multivariate normal distribution corresponding to the state variables. Let the changes in the scaled state variables x1, x2, x3, and x4 during the time interval (t, t + δt) be denoted by δx1, δx2, δx3, and δx4 respectively, where δxi(t) = xi(t + δt) − xi(t), 1, 2, 3, 4.

Under the assumptions of the original process on the sequence of transitions, we evaluate the mean vector and covariance matrix for δxi (i = 1, 2, 3, 4) during the time interval (t, t + δt) as follows. At first, we assume that we are in the state (S, I, C, A). So, the possible transitions from this state are:

(a) S increases by 1 at the rate μ
(b) S decreases by 1 at the rate μ
(c) S decreases by 1 and I increases by 1 at the rate βSI/N
(d) S decreases by 1 and I increases by 1 at the rate βηC SC/N
(e) S decreases by 1 and I increases by 1 at the rate βηA SA/N
(f) I decreases by 1 and C increases by 1 at the rate ϕI
(g) I decreases by 1 and A increases by 1 at the rate ρI
(h) I decreases by 1 at the rate μI
(i) C decreases by 1 and I increases by 1 at the rate ωC
(j) C decreases by 1 at the rate μC
(k) A decreases by 1 and I increases by 1 at the rate αA
(l) A decrease by 1 at the rate (μ + d)A.

The random variable δx1 equals $\frac{1}{N}$ in cases (a), $-\frac{1}{N}$ in cases (b), (c), (d), (e), and 0 in other cases. Similarly, δx2 equals $\frac{1}{N}$ in cases (c), (d), (e), (i), (k), $-\frac{1}{N}$ in cases (f), (g), (h), and 0 in other cases. δx3 equals $\frac{1}{N}$ in cases (f), $-\frac{1}{N}$ in cases (i), (j), and 0 in other cases. δx4 equals $\frac{1}{N}$ in cases (g), $-\frac{1}{N}$ in cases (k), (l), and 0 in other cases. Denoting x = (x1, x2, x3, x4)′ we have,

$\begin{array}{cc}& E\left(\delta \mathbit{x}\right)=b\left(\mathbit{x}\right)\delta t+o\left(\delta t\right)\end{array}$
$\begin{array}{c}\hfill \text{where}\phantom{\rule{4pt}{0ex}}b\left(\mathbit{x}\right)=\left(\begin{array}{c}\mu -\beta \left({x}_{2}+{\eta }_{C}{x}_{3}+{\eta }_{A}{x}_{4}\right){x}_{1}-\mu {x}_{1}\\ \beta \left({x}_{2}+{\eta }_{C}{x}_{3}+{\eta }_{A}{x}_{4}\right){x}_{1}-\left(\rho +\varphi +\mu \right){x}_{2}+\omega {x}_{3}+\alpha {x}_{4}\\ \varphi {x}_{2}-\left(\omega +\mu \right){x}_{3}\\ \rho {x}_{2}-\left(\alpha +\mu +d\right){x}_{4}\end{array}\right)\end{array}$

Now to derive the covariance matrix we need to find the Jacobian matrix at point x, which is obtained from b(x) as,

$\begin{array}{c}\hfill B\left(\mathbit{x}\right)=\frac{\partial b\left(\mathbit{x}\right)}{\partial \mathbit{x}}=\left(\begin{array}{cccc}-\beta \left({x}_{2}+{\eta }_{C}{x}_{3}+{\eta }_{A}{x}_{4}\right)-\mu & -\beta {x}_{1}& -\beta {\eta }_{C}{x}_{1}& -\beta {\eta }_{A}{x}_{1}\\ \beta \left({x}_{2}+{\eta }_{C}{x}_{3}+{\eta }_{A}{x}_{4}\right)& \beta {x}_{1}-\left(\rho +\varphi +\mu \right)& \beta {\eta }_{C}{x}_{1}+\omega & \beta {\eta }_{A}{x}_{1}+\alpha \\ 0& \varphi & -\left(\omega +\mu \right)& 0\\ 0& \rho & 0& -\left(\alpha +\mu +d\right)\end{array}\right)\end{array}$
Then we approximate B(x) at equilibrium point $\stackrel{^}{x}={\left({\stackrel{^}{\mathbit{x}}}_{1},{\stackrel{^}{x}}_{2},{\stackrel{^}{x}}_{3},{\stackrel{^}{x}}_{4}\right)}^{\prime }$ by $B\left(\stackrel{^}{\mathbit{x}}\right)$. So, covariance matrix of δx = (δx1, δ x2, δx3, δx4)′ is,

$V\left(\delta \mathbit{x}\right)=\frac{1}{N}S\left(\mathbit{x}\right)\delta t+o\left(\delta t\right)$ where,

$\begin{array}{c}\hfill S\left(\mathbit{x}\right)=\left(\begin{array}{cccc}\beta \left({x}_{2}+{\eta }_{C}{x}_{3}+{\eta }_{A}{x}_{4}\right){x}_{1}& -\beta {x}_{1}{x}_{2}& -\beta {\eta }_{C}{x}_{1}{x}_{3}& -\beta {\eta }_{A}{x}_{1}{x}_{4}\\ {\mu }^{\prime }+\mu {x}_{1}& & & \\ -\beta {x}_{1}{x}_{2}& \beta \left({x}_{2}+{\eta }_{C}{x}_{3}+{\eta }_{A}{x}_{4}\right){x}_{1}& -\varphi {x}_{2}-\omega {x}_{3}& -\rho {x}_{2}-\alpha {x}_{4}\\ & +\left(\rho +\varphi +\mu \right){x}_{2}+\omega {x}_{3}+\alpha {x}_{4}& & \\ -\beta {\eta }_{C}{x}_{1}{x}_{3}& -\varphi {x}_{2}-\omega {x}_{3}& \varphi {x}_{2}+\left(\omega +\mu \right){x}_{3}& 0\\ -\beta {\eta }_{A}{x}_{1}{x}_{4}& -\rho {x}_{2}-\alpha {x}_{4}& 0& \rho {x}_{2}\\ & & & +\left(\alpha +\mu +d\right){x}_{4}\end{array}\right)\end{array}$
We again approximate S(x) by $S\left(\stackrel{^}{\mathbit{x}}\right)$, where $\stackrel{^}{\mathbit{x}}$ is the equilibrium point. For large N, the process $\sqrt{N}\left(\mathbit{x}\left(t\right)-\stackrel{^}{\mathbit{x}}\right)$ is approximated by a multivariate Ornstein-Uhlenbeck (O-U) process [16], with a local drift matrix $B\left(\stackrel{^}{\mathbit{x}}\right)$ and local covariance matrix $S\left(\stackrel{^}{\mathbit{x}}\right)$.

The stationary distribution of this O-U process approximates the quasi-stationary distribution. It is approximately normal with mean zero and covariance matrix Σ, where Σ is obtained by solving

$\begin{array}{c}\hfill B\left(\stackrel{^}{\mathbit{x}}\right)\Sigma +\Sigma {B}^{\prime }\left(\stackrel{^}{\mathbit{x}}\right)=-S\left(\stackrel{^}{\mathbit{x}}\right).\end{array}$
Exact analytical solution for Σ is not straightforward [5]. However, since we are interested in calculating the CCS, we can easily solve the Eq (22) numerically using the parameter values and the equilibrium point.

Let σij be the solution of the (i, j)th element of Σ, where i, j = 1, …, 4. The diffusion approximation guides us to consider appropriate mean and variance-covariance matrix for the joint distribution of x(t). Thus we have,

$\begin{array}{c}\hfill \sqrt{N}\left(\mathbit{x}\left(t\right)-\stackrel{^}{\mathbit{x}}\right)\sim {N}_{4}\left(\mathbf{0},\Sigma \right),\phantom{\rule{4pt}{0ex}}\text{with}\phantom{\rule{4pt}{0ex}}\mathbit{x}\left(t\right)={\left({x}_{1}\left(t\right),{x}_{2}\left(t\right),{x}_{3}\left(t\right),{x}_{4}\left(t\right)\right)}^{\prime },\stackrel{^}{\mathbit{x}}={\left({\stackrel{^}{x}}_{1},{\stackrel{^}{x}}_{2},{\stackrel{^}{x}}_{3},{\stackrel{^}{x}}_{4}\right)}^{\prime },\end{array}$
$\begin{array}{c}\hfill \text{and}\phantom{\rule{4pt}{0ex}}\Sigma =\left(\begin{array}{cccc}{\sigma }_{11}& {\sigma }_{12}& {\sigma }_{13}& {\sigma }_{14}\\ {\sigma }_{21}& {\sigma }_{22}& {\sigma }_{23}& {\sigma }_{24}\\ {\sigma }_{31}& {\sigma }_{32}& {\sigma }_{33}& {\sigma }_{34}\\ {\sigma }_{41}& {\sigma }_{42}& {\sigma }_{43}& {\sigma }_{44}\end{array}\right)\end{array}$

The approximation for quasi-stationary distribution is obtained by using conditional truncated distributions of the above multivariate normal distribution. Thus in order to evaluate p•000, p•100, p•010, and p•001 we need to use a result from conditional truncated multivariate normal distribution, which is given below as theorem 1 (for proof see S1 Appendix). Even then, it is extremely difficult, if not impossible [5], to get an exact expression for q(d, μ). So, we evaluate an approximate expression of q(d, μ ) using Result 1 (for proof see S1 Appendix).

Theorem 1 Let YNp(μ, Σ) and write $\mathbit{y}={\left({\mathbit{y}}_{1}^{\prime },{\mathbit{y}}_{2}^{\prime }\right)}^{\prime }$, $\mathbit{\mu }={\left({\mathbit{\mu }}_{1}^{\prime },{\mathbit{\mu }}_{2}^{\prime }\right)}^{\prime }$, and $\Sigma =\left(\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right)$. Suppose instead of ${\mathbb{R}}^{p}$, Y is defined only on a truncated support c < y < d. Consider the partitions $\mathbit{c}={\left({\mathbit{c}}_{1}^{\prime },{\mathbit{c}}_{2}^{\prime }\right)}^{\prime }$ and $\mathbit{d}={\left({\mathbit{d}}_{1}^{\prime },{\mathbit{d}}_{2}^{\prime }\right)}^{\prime }$. Then, the conditional distribution of Y1 given y2 is given by

$\begin{array}{c}\hfill {f}^{*}\left({y}_{1}|{y}_{2}\right)=\frac{f\left({y}_{1}|{y}_{2}\right)}{{\int }_{{c}_{1}}^{{d}_{1}}f\left({y}_{1}|{y}_{2}\right)d{y}_{1}}\end{array}$
where f(y1|y2) is the conditional probability density function of Y1 given y2 i.e. ${\mathbit{Y}}_{1}|{\mathbit{y}}_{2}\sim {N}_{q}\left({\mathbit{\mu }}_{1}+{\Sigma }_{12}{\Sigma }_{22}^{-1}\left({\mathbit{y}}_{2}-{\mathbit{\mu }}_{2}\right),{\Sigma }_{11.2}={\Sigma }_{11}-{\Sigma }_{12}{\Sigma }_{22}^{-1}{\Sigma }_{21}\right)$

Result 1 Let $\varphi \left(x\right)=\frac{1}{\sqrt{2\pi }}{e}^{-{x}^{2}/2}$ and $\Phi \left(x\right)={\int }_{-\infty }^{x}\varphi \left(t\right)dt\phantom{\rule{4pt}{0ex}}for\phantom{\rule{4pt}{0ex}}any\phantom{\rule{4pt}{0ex}}\mathbit{x}\in \left(-\infty ,\infty \right)$. Then, an approximate expression of ${q}_{•}^{\left(d,\mu \right)}$ is given as,

$\begin{array}{c}\hfill {q}_{•}^{\left(d,\mu \right)}=\frac{\mu {p}_{•100}+\mu {p}_{•010}+\left(d+\mu \right){p}_{•001}}{1-{p}_{•000}}\end{array}$
$\begin{array}{cc}\hfill where& {p}_{•100}\approx \frac{1}{2N\sqrt{{\sigma }_{22}^{*}}}\frac{\varphi \left(\frac{{\stackrel{^}{x}}_{2}^{*}}{\sqrt{{\sigma }_{22}^{*}}}\right)}{\Phi \left(\frac{{\stackrel{^}{x}}_{2}^{*}}{\sqrt{{\sigma }_{22}^{*}}}\right)}.\frac{1}{2N\sqrt{{\sigma }_{33}^{*}}}\frac{\varphi \left(\frac{{\stackrel{^}{x}}_{3}^{*}}{\sqrt{{\sigma }_{33}^{*}}}\right)}{\Phi \left(\frac{{\stackrel{^}{x}}_{3}^{*}}{\sqrt{{\sigma }_{33}^{*}}}\right)}.\frac{1}{2N\sqrt{{\sigma }_{44}}}\frac{\varphi \left(\frac{{\stackrel{^}{x}}_{4}}{\sqrt{{\sigma }_{44}}}\right)}{\Phi \left(\frac{{\stackrel{^}{x}}_{4}}{\sqrt{{\sigma }_{44}}}\right)}\\ & {p}_{•010}\approx \frac{1}{2N\sqrt{{\sigma }_{33}^{**}}}\frac{\varphi \left(\frac{{\stackrel{^}{x}}_{3}^{**}}{\sqrt{{\sigma }_{3}^{**}}}\right)}{\Phi \left(\frac{{\stackrel{^}{x}}_{3}^{**}}{\sqrt{{\sigma }_{33}^{**}}}\right)}.\frac{1}{2N\sqrt{{\sigma }_{22}^{**}}}\frac{\varphi \left(\frac{{\stackrel{^}{x}}_{2}^{**}}{\sqrt{{\sigma }_{2}^{**}}}\right)}{\Phi \left(\frac{{\stackrel{^}{x}}_{2}^{**}}{\sqrt{{\sigma }_{22}^{**}}}\right)}.\frac{1}{2N\sqrt{{\sigma }_{44}}}\frac{\varphi \left(\frac{{\stackrel{^}{x}}_{4}}{\sqrt{{\sigma }_{44}}}\right)}{\Phi \left(\frac{{\stackrel{^}{x}}_{4}}{\sqrt{{\sigma }_{44}}}\right)}\\ & {p}_{•001}\approx \frac{1}{2N\sqrt{{\sigma }_{44}^{***}}}\frac{\varphi \left(\frac{{\stackrel{^}{x}}_{4}^{***}}{\sqrt{{\sigma }_{44}^{***}}}\right)}{\Phi \left(\frac{{\stackrel{^}{x}}_{4}^{***}}{\sqrt{{\sigma }_{44}^{***}}}\right)}.\frac{1}{2N\sqrt{{\sigma }_{22}^{***}}}\frac{\varphi \left(\frac{{\stackrel{^}{x}}_{2}^{***}}{\sqrt{{\sigma }_{22}^{***}}}\right)}{\Phi \left(\frac{{\stackrel{^}{x}}_{2}^{***}}{\sqrt{{\sigma }_{22}^{***}}}\right)}.\frac{1}{2N\sqrt{{\sigma }_{33}}}\frac{\varphi \left(\frac{{\stackrel{^}{x}}_{3}}{\sqrt{{\sigma }_{33}}}\right)}{\Phi \left(\frac{{\stackrel{^}{x}}_{3}}{\sqrt{{\sigma }_{33}}}\right)}\end{array}$
where ${\stackrel{^}{x}}_{i}^{*},{\sigma }_{ii}^{*}\phantom{\rule{4pt}{0ex}}for\phantom{\rule{4pt}{0ex}}i=2,3$, ${\stackrel{^}{x}}_{i}^{**},{\sigma }_{ii}^{**}\phantom{\rule{4pt}{0ex}}for\phantom{\rule{4pt}{0ex}}i=2,3$, ${\stackrel{^}{x}}_{i}^{***},{\sigma }_{ii}^{***}\phantom{\rule{4pt}{0ex}}for\phantom{\rule{4pt}{0ex}}i=2,4$ are obtained from the truncated conditional distribution of multivariate normal distribution as given in (23).

Once we find ${q}_{•}^{\left(d,\mu \right)}$, we have an expression for expected time to extinction $\stackrel{^}{E}\left({\tau }_{Q}\right)$ using Eq (12). However, this involves N, which is unknown. Now, our aim is to find the above expression in another way so that we can equate them and solve for N . So, following [6], if we can find the quasi-period (${\stackrel{^}{T}}_{0}$), we may use the relation $\stackrel{^}{E}\left({\tau }_{Q}\right)\phantom{\rule{4pt}{0ex}}\text{log}\phantom{\rule{4pt}{0ex}}2={\stackrel{^}{T}}_{0}$ to get the value of N, which is the critical community size.

We determine the quasi-period of the oscillation about the critical point using linearisation method [17]. Note that for our model, the linearised system about the equilibrium point $\stackrel{^}{\mathbit{x}}={\left({\stackrel{^}{x}}_{1},{\stackrel{^}{x}}_{2},{\stackrel{^}{x}}_{3},{\stackrel{^}{x}}_{4}\right)}^{\prime }$ can be written as:

$\begin{array}{c}\hfill \frac{d{\mathbit{x}}^{*}}{dt}=\left(\begin{array}{cccc}-\beta \left(1+{\eta }_{C}\frac{\varphi }{{\xi }_{2}}+{\eta }_{A}\frac{\rho }{{\xi }_{1}}\right){\stackrel{^}{x}}_{2}-\mu & -\frac{\beta }{{R}_{0}}& -\frac{\beta {\eta }_{C}}{{R}_{0}}& -\frac{\beta {\eta }_{A}}{{R}_{0}}\\ \beta \left(1+{\eta }_{C}\frac{\varphi }{{\xi }_{2}}+{\eta }_{A}\frac{\rho }{{\xi }_{1}}\right){\stackrel{^}{x}}_{2}& \frac{\beta }{{R}_{0}}-\left(\rho +\varphi +\mu \right)& \frac{\beta {\eta }_{C}}{{R}_{0}}+\omega & \frac{\beta {\eta }_{A}}{{R}_{0}}+\alpha \\ 0& \varphi & -\left(\omega +\mu \right)& 0\\ 0& \rho & 0& -\left(\alpha +\mu +\alpha \right)\end{array}\right){\mathbit{x}}^{*}\end{array}$
where ${\mathbit{x}}^{*}=\mathbit{x}-\stackrel{^}{\mathbit{x}}$.

Now we can find the eigenvalues of the matrix in Eq (27) and hence the angular frequency, provided there are imaginary roots. Using Eq (27) and putting the values of the parameters, we can also find the angular frequency (θ) numerically. Hence the quasi-period (${\stackrel{^}{T}}_{0}$, say) is obtained as ${\stackrel{^}{T}}_{0}=\frac{2\pi }{\theta }$, which is independent of N. Now noting the relation $\stackrel{^}{E}\left({\tau }_{Q}\right)\phantom{\rule{4pt}{0ex}}\text{log}\phantom{\rule{4pt}{0ex}}2={\stackrel{^}{T}}_{0}$ and using Result 1, we can solve for N.

Since we are dealing with a system consisting of four equations, it is not possible to obtain an explicit or compact expression for CCS (i.e. N) in terms of the model parameters. Hence we find an approximate value of CCS numerically.

Our approach is not restricted to only two [6] or three variables [5]. It generalizes the numerical calculation of CCS for any finite number of variables. Thus, our method for finding CCS is a general one and is applicable to any number of variables. Let n be the number of variables and m be the number of equations in the system of differential equations that explains the disease dynamics. An algorithm to calculate CCS for such a process is given below.

• Step I: Calculate basic reproduction number and endemic equilibrium for the system of differential equations i.e. the model.
• Step II: Write down a fully stochastic model and the joint distribution of the state variables conditioning on non-extinction of the disease.
• Step III: Write the expression for expected time to extinction based on quasi-stationary distribution.
• Step IV: Introduce diffusion approximation to obtain the local drift matrix and local covariance matrix from multivariate Ornstein-Uhlenbeck process.
• Step V: Using conditional truncated multivariate normal distribution, calculate expected time to extinction ($\stackrel{^}{E}\left({\tau }_{Q}\right)$, say)
• Step VI: Using linearisation of the model, calculate the angular frequency (${\stackrel{^}{T}}_{0}$, say).
• Step VII: Merge expected time to extinction with angular frequency by the relation $\stackrel{^}{E}\left({\tau }_{Q}\right)\phantom{\rule{2pt}{0ex}}\text{log}\phantom{\rule{2pt}{0ex}}2={\stackrel{^}{T}}_{0}$, and solve for N, which is the CCS with respect to the model.

Note that CCS thus calculated is an approximate one as we have made some approximations at few stages. However, this algorithm is a general one and may be applied to any model involving a system of differential equations that explains the disease dynamics. However, the basic reproduction number (R0) must be greater than 1 in order to calculate CCS using our method.

## Results

### Case study

For HIV transmission, we can calculate the critical community size for any region or community. Note that this value of CCS and time to extinction are approximate as we have made some approximation while applying diffusion approximation using conditional truncated multivariate normal distribution. All calculations are based on the assumed model (1)-(4). We suppose that an individual in the susceptible compartment/class (S) either moves to the infected (I) compartment or infected under treatment (C) compartment or infected with AIDS (A) compartment; while an individual from infected (I) class may receive treatment and move to infected under treatment (C) class or get infected with AIDS and move to infected with AIDS (A) class; while an individual after treatment may have reduced viral load and thus move from class C to the infected (I) class; similarly, an individual infected with AIDS (A) may have reduced viral load and move to infected (I) class. Death may occur in each class while immigration occurs only in the susceptible class.

It is clear that in order to calculate CCS for a region or community, we need to know the values of the parameters in the model. These parameter values are usually estimated based on different studies. As an illustration, we consider the values of the parameters based on the data available for Madagascar. However, the values of all parameters are not available in the literature. In such cases, we have assumed a few values considering the meaning and range of the parameters. Based on the available data for 2018 from the UNAIDS (www.unaids.org/en/regionscountries/countries/madagascar), we see that the number of people living with HIV who are on ART is 3500 whereas the number of adults and children living with HIV is 39000. The number of adults and children newly infected with HIV is 6100 and the number of deaths due to AIDS is 1700. Based on this information and following Nåsell [6] we have considered β = 0.156 × 365 = 56.94 and d = 0.044 × 365 = 16.06. Based on the available reports (http://cfs.hivci.org/country-factsheet.html) from the World Health Organisation (WHO), we have ω = 0.03 × 365 = 10.95, ϕ = 0.0035 × 365 = 1.2775, ηC = 0.000105, and μ = 0.015. Other parameters such as α = 0.001 × 365 = 0.365, ηA = 1.1, and ρ = 0.98 × 365 = 357.7 are based on available literature and intuitive meaning of the parameters. Given the poor socio-economic conditions of Madagascar, the 90-90-90 HIV treatment target initiative by UNAIDS has not performed remarkably well. At severe stages of HIV infection, available medication has very little effect. This information has led us to opt for such a high value of ρ and small value for α. Here, β is the contact rate (per year) of an HIV-infected person with the susceptible pool of individuals and ρ signifies the rate of progress (per year) of HIV-infected individuals to become AIDS patient. Similarly, we define other mentioned rates.

For this set of parameters, the CCS is obtained as 4585. So we can say that if the population in a commune drops below 4585, the infection will die out automatically. Obviously, the range of CCS for the HIV epidemic is far less than comparatively more infectious diseases such as measles that has much higher incidence rate. UNICEF reports that there were 244,607 cases of measles and 1,080 deaths in Madagascar from August 2018 to November 2019 (www.unicef.org/press-releases/measles-outbreaks-continue-unabated-five-countries-accounted-nearly-half-all-measles). Naturally, Nåsell [6] obtained CCS as 561, 000 (with R0 = 14) and Bartlett [4] estimated CCS as 250, 000 to 300, 000 in England and Wales for measles. Note that, using the above parameters for Madagascar and assuming homogenous mixing within communities, we obtain R0 as 4.059, which essentially means that in a completely susceptible population, a single infection may produce 4.059 secondary infections in the duration of infectiveness. Our calculation also indicates that temporary eradication or fade-out could happen after 7.365 years unless the infection is reintroduced from outside. This fade-out tends to happen in local spatial regions until the reintroduction of the disease occurs from other areas. However, this does not imply that a person already infected with HIV will be removed from the population within this period. In fact, she/he may continue to exist for more than 7.365 years and we assume that no new infection will be spread by her/him.

Next, we study the nature of CCS by varying one parameter and keeping others fixed at the above values. For this we consider the most important parameters (β, d, ϕ, ω, α, and ρ ) in our model. Fig 1 indicates the relation of CCS with different parameters. Here, β is the contact rate for HIV transmission. An increase of β signifies an increment of infection circulation in the population. So the susceptible population will get infected at a much higher rate. In order to stop the infection, the population size should be very small. So if the number of the susceptible population falls below some threshold, the chance of the spread of infection will diminish. This indicates that smaller the population size, lesser would be the chance of spreading the infection. Thus, it is expected that if β increases, CCS will decrease (Fig 1A). d is the rate of death from class A. So as d increases, the spread of infection will decrease. For large values of d, a large susceptible population will remain unaffected. Hence the CCS value will be high for large d (Fig 1B). ϕ and ω respectively are the rates at which individuals in class I move to C and individuals in class C move to I. So, as ϕ increases, the number of individuals receiving treatment increases. It is expected that individuals under treatment will spread less infection than those who are infected but not treated. So, as more individuals go to class I from class C, the chance of spreading infection in the susceptible population decreases. But, it is known that the infection stays in a community long enough to engulf the entire susceptible population. Therefore, the infection might be eradicated if the susceptible population decreases. This would mean that, under the condition that ϕ increases, eradication of infection is expected to happen as the CCS value decreases (Fig 1C). On the contrary, as ω increases, CCS also increases (Fig 1D). Two other parameters viz. α and ρ represent the rates at which individuals in class A move to I and individuals in class I move to A respectively. For similar reasons, as α increases, CCS decreases (Fig 1E) but as ρ increases, CCS increases (Fig 1F).

Fig 1
Relation of CCS with different parameters.

## Discussion

Bartlett [3, 4] describes fade-out tendency of measles epidemic for small populations. He noticed that for such populations as the size of the susceptible pool of individuals is equal to CCS, measles has a 50% chance of fading out. Bartlett observed that for some cities in the United States where the population size is below CCS, the epidemic ceases to exist in absence of reintroduction of the infection from outside. But for larger cities, the tendency of damping down of successive epidemics is offset by random variability. The idea of epidemic fade-out in smaller spatial locations provides a clue to systematically treat epidemic rather than concentrating treatment for an entire continent or a large community at once.

Shocking statistics (www.prb.org/thestatusofthehivaidsepidemicinsubsaharanafrica) reveal that 16 countries in Sub-Saharan Africa constitute nearly 4% of the world’s population but account for more than 50% of HIV infections worldwide. In some of these countries, the prevalence rate is up to 30%. However, efforts to control HIV epidemic situations have resulted in a decrease in the disease prevalence in some African countries. Total eradication of highly prevalent diseases such as HIV infection in the context of countries in Sub-Saharan Africa is almost impossible in a very short time period. So we may instill the idea of CCS to systematically eradicate HIV. If we concentrate on the treatment of HIV in a more focussed manner in smaller spatial regions, it may accelerate the HIV annihilation process unless the disease is reintroduced from outside.

When we look at the geographical locations of Uganda, Malawi, Eswatini, and Madagascar, interestingly we find that being at the heart of the Great Lakes region, Uganda is surrounded by three lakes among which one is a fresh water lake. These lakes, to some extent separate Uganda from the rest of the African countries. Malawi again is separated by a freshwater lake (Lake Malawi) that comprises of 25% of its area. But Eswatini is a landlocked country in Southern Africa. We observe that, even with quite impressive progress towards achieving the 90-90-90 target, countries that are more landlocked have higher rates of incidence and prevalence compared to other countries that are isolated to varying extent. Although being in the AIDS belt of Africa, Madagascar has shown remarkably lower rates of incidence and prevalence even with little success of the 90-90-90 target. It is possibly because of the isolated location of Madagascar that separates it from the mainland of the African continent to a great extent.

Metcalf et al. [8] pointed out with respect to some childhood infections that the relative risk of extinction of infections on islands is almost double that in the mainland. Such infections cease to exist at a significantly higher rate in islands where access to/from the mainland is much less. But the persistence of infection within a region in the mainland is fuelled constantly by recolonization. Hence, although the epidemic infection should apparently die out because of the natural stochastic extinction property (damping down) of itself after invading the total susceptible subpopulation in the mainland, recolonization introduces infection and hinders the damping down process. Thus, our observations are consistent with the views of this paper.

Madagascar has 111 large and small cities or districts with 39.6% of them have population size less than 20000, 30.6% with population size somewhere between 20000 − 30000, 16.2% has 30000 − 40000 and 5.4% have population more than 100000 (http://worldpopulationreview.com/countries/madagascar-population/cities/). Further, each district is comprised of several communes of an even smaller population. So our idea is that if the population size falls below the CCS of Madagascar which is 4585 in each spatially separated population, HIV infection may cease to exist in different spatial locations, after a mean time of 7.365 years, provided infection is not re-introduced from outside. But this will definitely not happen without continuing the controlling measures like educating people about the disease and taking necessary preventive measures/treatment available. Through this paper, we are emphasizing that the systemic spreading of control measures could lead to an acceleration in achieving the 90-90-90 goal of UNAIDS.

We observe from our case study that in addition to the existing controlling measures of HIV epidemics in Madagascar, employing the idea of CCS might accelerate the eradication of the infection from smaller populations after the corresponding mean extinction time. Besides, CCS could be looked upon as an index of performance of epidemic controlling strategies in different countries.

## Conclusions

Combating pandemic and epidemic outbreaks require efficient controlling strategies, as the available treatment facilities and/or vaccines become limited in such situations. But, the infection persists in the community long enough to engulf the entire susceptible population. Hence, in large communities, the infection recurs in the subsequent timeline. Interestingly, Bartlett observed that in smaller communities the infection dies out if the susceptible population drops below a threshold size (CCS). Later, Nåsell formulated CCS for infectious diseases having Susceptible-Infected model dynamics. But the analytical calculation of CCS for more complex disease dynamics or higher-order models are quite cumbersome. So, we propose a simplified computation based approach for rapid calculation of the approximate value of CCS, for HIV disease dynamics with four differential equations. Since our method is free from any restriction on the number of variables in the system, it could be applied to other diseases with higher-order model dynamics. We envisage CCS for devising control strategies during epidemic outbreaks such as COVID-19 [18]. Because, during epidemic outbreaks with limited treatment facilities, keeping the susceptible population below CCS by enforcing appropriate containment plans might control the disease substantially.

So to sum up, through this paper, firstly, we have generalized the approach for CCS calculation for a system of differential equations having more than two variables. Our method for a system containing four variables may be extended to higher dimensions for any infectious disease dynamics. This is because our method is able to tackle multivariate normal distribution for calculating expected time to extinction as a function of CCS (Result 1). Secondly, we have developed a method to unleash the potential that CCS holds in devising control strategies by promptly determining its approximate numerical value from any disease dynamics using a simplified computation approach. In particular, we believe that our method will provide a suitable measure for accelerated extermination of epidemic infections by aiding policymakers in making decisions on control strategies. Thirdly, we hypothesize that, if we start treating spatially isolated regions or smaller communes and then progress to the district level, and after that to larger regions systematically, the annihilation of HIV infection could be accelerated. Although it is way difficult to reach out to smaller communes with treatments and other controlling measures (such as educating them about the infection etc.), CCS could be worth giving a trial, as this could utilize the natural force of damping down of epidemic oscillations to eradicate HIV infection from smaller populations.

## References

SoperHE. The interpretation of periodicity in disease prevalence. Journal of the Royal Statistical Society. 1929;92(1):3473. doi: 10.2307/2341437

Bartlett M. S. 1956 Deterministic and stochastic models for recurrent epidemics. In: Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability (ed. J. Neyman); 1956. p. 81–109.

BartlettMS. Measles periodicity and community size. Journal of the Royal Statistical Society Series A (General). 1957;120(1):4870. doi: 10.2307/2342553

BartlettM. The critical community size for measles in the United States. Journal of the Royal Statistical Society: Series A (General). 1960;123(1):3744. doi: 10.2307/2343186

AnderssonH, BrittonT. Stochastic epidemics in dynamic populations: quasi-stationarity and extinction. Journal of mathematical biology. 2000;41(6):559580.

NåsellI. A new look at the critical community size for childhood infections. Theoretical population biology. 2005;67(3):203216.

NåsellI. Stochastic models of some endemic infections. Mathematical biosciences. 2002;179(1):119.

MetcalfCJE, HampsonK, TatemAJ, GrenfellBT, BjørnstadON. Persistence in epidemic metapopulations: quantifying the rescue effects for measles, mumps, rubella and whooping cough. PloS one. 2013;8(9):e74696

SilvaCJ, TorresDF. Global stability for a HIV/AIDS model. Communications Faculty of Sciences University of Ankara Series A1 Mathematics and Statistics. 2018;67:93101.

10

SilvaCJ, TorresDF. A TB-HIV/AIDS coinfection model and optimal control treatment. Discrete & Continuous Dynamical Systems A. 2015;35(9):46394663.

11

DjordjevicJ, SilvaCJ, TorresDF. A stochastic SICA epidemic model for HIV transmission. Applied Mathematics Letters. 2018;84:168175.

12

ArifMS, RazaA, RafiqM, BibiM. A reliable numerical analysis for stochastic hepatitis B virus epidemic model with the migration effect. Iranian Journal of Science and Technology, Transactions A: Science. 2019;43(5):24772492.

13

ArifMS, RazaA, RafiqM, BibiM, FayyazR, NazM, et al A reliable stochastic numerical analysis for typhoid fever incorporating with protection against infection. Computers, Materials & Continua. 2019;59(3):787804.

14

RazaA, ArifMS, RafiqM. A reliable numerical analysis for stochastic dengue epidemic model with incubation period of virus. Advances in Difference Equations. 2019;2019(1):32

15

HeffernanJM, SmithRJ, WahlLM. Perspectives on the basic reproductive ratio. Journal of the Royal Society Interface. 2005;2(4):281293.

16

NåsellI. On the time to extinction in recurrent epidemics. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 1999;61(2):309330.

17

DietzK. Transmission and control of arbovirus diseases. Epidemiology. 1975;104:121.

18

DasS, GhoshP, SenB, PyneS, MukhopadhyayI. Critical community size for COVID-19: A model based approach for strategic lockdown policy. Statistics and Applications. 2020;18(1):181196.
This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Citing articles via
https://www.researchpad.co/tools/openurl?pubtype=article&doi=10.1371/journal.pone.0244543&title=Determination of critical community size from an HIV/AIDS model&author=&keyword=&subject=Research Article,Medicine and health sciences,Epidemiology,HIV epidemiology,Medicine and health sciences,Medical conditions,Infectious diseases,Viral diseases,HIV infections,People and Places,Geographical Locations,Africa,Madagascar,Medicine and Health Sciences,Epidemiology,Infectious Disease Epidemiology,Medicine and Health Sciences,Medical Conditions,Infectious Diseases,Infectious Disease Epidemiology,Medicine and Health Sciences,Epidemiology,Disease Dynamics,Medicine and Health Sciences,Medical Conditions,Infectious Diseases,Infectious Disease Control,Biology and Life Sciences,Microbiology,Medical Microbiology,Microbial Pathogens,Viral Pathogens,Immunodeficiency Viruses,HIV,Medicine and Health Sciences,Pathology and Laboratory Medicine,Pathogens,Microbial Pathogens,Viral Pathogens,Immunodeficiency Viruses,HIV,Biology and Life Sciences,Organisms,Viruses,Viral Pathogens,Immunodeficiency Viruses,HIV,Biology and Life Sciences,Organisms,Viruses,Immunodeficiency Viruses,HIV,Biology and life sciences,Organisms,Viruses,RNA viruses,Retroviruses,Lentivirus,HIV,Biology and Life Sciences,Microbiology,Medical Microbiology,Microbial Pathogens,Viral Pathogens,Retroviruses,Lentivirus,HIV,Medicine and Health Sciences,Pathology and Laboratory Medicine,Pathogens,Microbial Pathogens,Viral Pathogens,Retroviruses,Lentivirus,HIV,Biology and Life Sciences,Organisms,Viruses,Viral Pathogens,Retroviruses,Lentivirus,HIV,Physical Sciences,Mathematics,Probability Theory,Probability Distribution,Normal Distribution,