Modelling the effects of the contaminated environments on tuberculosis in Jiangsu, China.

Tuberculosis (TB) is still an important public health issue in Jiangsu province, China. In this study, based on the TB transmission routes and the statistical data of TB cases, we formulate a novel TB epidemic model accounting for the effects of the contaminated environments on TB transmission dynamics. The value of this study lies in two aspects. Mathematically, we define the basic reproduction number, R0, and prove that R0 can be used to govern the threshold dynamics of the model. Epidemiologically, we find that the annual average R0 is 1.13,>1 and TB in Jiangsu is an endemic disease. Therefore, in order to control the TB in Jiangsu efficiently, we must decrease the virus shedding rate or increase the recovery rates, and increase the environmental clearance rate.

Introduction

Tuberculosis (TB) is an ancient and chronic infectious disease, which is caused by infection with the Mycobacterium tuberculosis (MTB) (Blower et al., 1996). Normally, the TB bacteria are put into the air when a person with TB disease, and TB patients are mainly transmitted by droplets produced by coughing, sneezing, laughing, loud talking, etc. Droplet transmission is the most important way of transmission of TB. When a person breathes in TB bacteria, the bacteria can settle in the lungs and begin to grow (Global, 2018).

It is now widely believed that droplet transmission occurs when a person is in close contact (within 1 m) with someone who has respiratory symptoms (e.g., coughing or sneezing) and is therefore at risk of having his/her mucosae (mouth and nose) or conjunctiva (eyes) exposed to potentially infective respiratory droplets. Transmission may also occur through fomites in the immediate environment around the infected person (Ong et al., 2020). MTB is so small that normal air currents can keep the particles containing MTB airborne and transport them through rooms or some buildings (Dye and Williams, 2000). Thank the insightful work on HFMD (Wang et al., 2016a, Wang et al., 2016b) and COVID-19 (Ding et al., 2020, World Health, 2020), we can believe that MTB can attach to things (such as door handles, towels, handkerchiefs, toys, utensils, bed and toilet seat, bathroom washbasin tap lever, bathroom ceiling-exhaust louvre and stethoscope or thermometer, and so on) used by the TB patients.

TB is closely associated with overcroding and malnutrition, which makes it to be one of the major diseases in poor areas (Lawn and Zumla, 2011). Those at high risk thus include: people who inject illicit drugs, inhabitants and employees of locales where vulnerable people gather (e.g. prisons and homeless shelters), medically underprivileged and resource-poor communities, high-risk ethnic minorities, children in close contact with high-risk category patients, and health-care providers serving these patients (Griffith and Kerr, 1996), include alcoholism (Lawn and Zumla, 2011) and diabetes mellitus (threefold increase) (Restrepo, 2007).

Currently, there are approximately 95 of the estimated 8 milion new cases of TB occuring in developing countries each year, and two-thirds of which appears in India (27%), China (9%), Indonesia (8%), Philippines (6%), Pakistan (5%), Nigeria (4%), Bangladesh (4%) and South Africa (3%), where 80 occur among people between the ages of 15 to 59 years. And only 6% of global cases occur in Europe (3%) and the Americas (3%) (Global, 2018).

In recent years, the Chinese government has increased its investment in public health, and the laws and regulations for disease prevention and control has been constantly improved, which provides an important basic guarantee for coping with major public health emergencies, preventing infectious diseases and ensuring the health of the people, thus effectively controlling major infectious diseases. In particular, the incidence of tuberculosis has decreased significantly. According to the report from the World Health Organization, in 2011, China had an estimation of 1.4 million existing TB and 1 million incident TB; in 2017, China had an estimation of 778,390 existing TB and 773,150 incident TB (Global, 2018); in 2018, China had 823,342 new and relapse TB (Survey of the epidemic situation of notifiable infectious diseases in China, 2018). Obviously, TB is still an important public health issue in China (Hu and Sun, 2013).

It was known that 80% of TB exists in rural areas, particularly in north and north-western regions with low socioeconomic status in China (Hu and Sun, 2013). But in the past 20 years in Jiangsu province, China (see Fig. 1 for the location of Jiangsu province in China), one of the most developed areas in China in economy, technology and culture and the total output is one of the largest in the nation, TB ranked first in the number of notifiable B infectious diseases (The reported tuberculosis cases in Jiangsu province, 2018). In 2011, Jiangsu had 39,589 existing TB, and in 2017, 28,402 existing TB and in 2018, there is 26,506 incident TB. In particular, the incidence of tuberculosis has decreased significantly (see Fig. 3(a) for more details). However, the situation of prevention and control of TB in Jiangsu province is still very serious.

Fig. 1

The location of Jiangsu province in China.

Fig. 3

The TB incidences times series and the relative infectivity in Jiangsu, China from 2009 to 2018. Panel (a) shows the monthly number of TB incidences. Panel (c) is an annualised version of panel (a). Panel (b) shows the relative TB infectivity from 2009 to 2018. Panel (d) is an annualised version of panel (b), where the short bars are the relative infectivity of each month of different years, and the diamonds are the average of each month.

It is worthy to notice that mathematical models have played a key role in the formulation of TB control strategies. Waaler et al. (1962) introduced the first mathematical model for TB in ordinary differential equations. The simplest TB transmission models include classes of susceptible, exposed, and infectious individuals, and hence, they are known as the SEI models. Of course, there are more factors that includes drug-resistant strains, fast and slow progression, confection with HIV, relapse, reinfection, migration, treatment, seasonality, and vaccination are incorporating into studying the transmission dynamics has been searched by many authors. Dye et al. (1998) present a model with explicit fast and slow progression from two latent classes. Ziv et al. (2004) used mathematical models to predict the potential public health impact of new TB vaccines in high-incidence countries. Porco and Blower (1998) included the aspect of disease relapse into their model. Then, there are mathematical models of TB including reinfection in some authors’ article and they assumed that the rate of reinfection is a multiple of the rate of first infection (Feng et al., 2000, Suzanne et al., 2005, Liu et al., 2010).

On the other hand, the TB can survive for a long period outside the host in suitable conditions, and hence contaminated environments may play important role in TB infection. There are some scholars investigated the disease dynamics in the contaminated environments. Wang et al., 2016a, Wang et al., 2016b, Chadsuthi and Wichapeng (2018) investigated the roles that asymptomatic individuals and contaminated environments played in Hand-foot-mouth disease dynamics. Machado et al. (2017) developed and implemented an integrative epidemiologic cross-sectional study that allows identifying and characterising exposure pathways of populations living and working on the shores of a contaminated estuarine environment.

There naturally comes a question that how do the contaminated environments affect the transmission dynamics of TB in Jiangsu, China?.

The main focus of this paper is to investigate how contaminated environments affect TB dynamics through studying the threshold dynamics of a general TB model. And the rest of the paper is organized as follows. In Section 2, we formulate the model in details. In Section 3, we give the dynamics analysis of the model, we introduce the basic reproduction number and prove that can be used to govern the threshold dynamics of the model. In Sections 4, we give the TB epidemics in Jiangsu, China via numerical simulations. In Section 5, we provide a brief discussion and the summary of the main results.

Model derivation

Suppose that the total population individuals divide into susceptible , infectious but not yet symptomatic, i.e., pre-symptomatic , infectious with symptoms , and recovered . We further consider the TB virus concentration in environment as , which is the density of pathogen of the contaminated environments including door handles, towels, handkerchiefs, toys, utensils, bed and toilet seat, bathroom washbasin tap lever, bathroom ceiling-exhaust louvre etc. at time t. And our model involves two typical transmissions: one is the direct transmission between susceptible and infected individuals (including pre-symptomatic and symptomatic ) with rate of ; the other is the indirect transmission to susceptible individuals and infected individuals (i.e., and ) by contaminated environments with rate of . A seasonality in the long-term patterns of TB incidences time series can be observed evidently (see Fig. 3), and there is a growing awareness that seasonality can cause population fluctuations ranging from annual cycles to multiyear oscillations (Liu et al., 2010), and hence we assume that the transmission rates and to be continuous and non-negative periodic functions with period of . A flow diagram describing the model is depicted in Fig. 2 .

Fig. 2

Flow diagram representing TB transmission routes.

Thus we can establish the following TB epidemic model involving five ordinary differential equations:with the initial conditionsand .The meanings of each variables and parameters in model (2.1) are as follows. with initial conditions

: the recruitment rate of susceptible;

: the per capita natural mortality rate;

: the rate of the susceptible get infected by direct individual-to-individual transmission;

: the rate of the susceptible get infected by indirect contaminated environment transmission;

p: proportion of TB symptomatic infectious by direct transmission;

q: proportion of TB symptomatic infectious by indirect transmission;

: reactivation rate the pre-symptomatic infectious;

: the recovery rate of the pre-symptomatic infectious;

: the recovery rate of the symptomatic infectious;

: the rate from recovered to susceptible.

: disease-related death;

: the virus shedding rate from symptomatic infected individuals;

c: the clearance rate of the virus in the environments.For the sake of conveniently analysis, set . Then model (2.1) with (2.2) is equivalent to the following system:

Model (2.3) has a unique and bounded solution with the initial value , i.e,

Considering the non-negativity of I, i.e., , the total population can be determined by the following model:

It is easy to see that the linear differential equation has a unique equilibrium , which is globally asymptotically stable. The comparison principle implies that .Then , there exists a such that

Then we haveIt follows from the comparison principle that . This completes the proof.

Let be the fundamental solution matrix of equation , where is a continuous, cooperative, irreducible and -periodic functional matrix. Let be the spectral radius of .

[Zhang and Zhao, 2007, Lemma 2.1] Let . Then there exists a positive -periodic function such that is a solution of .

Dynamics analysis

Basic reproduction number

Following (Diekmann et al., 1990, Van den Driessche and Watmough, 2002, Diekmann and Heesterbeek, 2000), let be the input rate of newly infected individuals and be the rate of transfer of individuals, thenand

Obviously model (2.3) admits a disease free equilibrium (DFE) .

Then

Let be the evolution operator of the linear -periodic system

That is, for each , the matrix satisfieswhere is the identity matrix. Thus, the monodromy matrix of (3.1) equals .

Following the method established by Wang and Zhao (2008), let be -periodic in s and the initial distribution of infectious individuals. So is the rate of new infections produced by the infected individuals who are introduced at time s. When gives the distribution of those infected individuals who are newly infected by and remain in the infected compartments at time t. Naturally,is the distribution of accumulative new infections at time t produced by all those infected individuals introduced at time previous to t.

Let be the ordered Banach space of all -periodic functions from to , which is equipped with the maximum norm and the positive cone . Then we can define a linear operator implies thatwhich is called the next infection operator, and the spectral radius of is defined as the basic reproduction number:

In order to characterise , we introduce the linear -periodic systemwith parameter . Let be the evolution operator of system (3.3) on . Clearly, . Hence, we derive

Following the general calculation procedure in Wang and Zhao (2008, Theorem 2.1), the basic reproduction number is the unique solution of .

Following Wang and Zhao (2008), we can obtain the relation between and shown in the following lemma.

[Wang and Zhao, 2008, Theorem 2.2] The following statements are valid:

(i)  iff ;

(ii)  iff ;

(iii)  iff .

In the special case of and , the basic reproduction number of model (2.3) is , i.e.,whereHere, indicates the average number of secondary infections generated by a single infected individual introduced into a completely susceptible population directly during their life cycle. indicates the average number of secondary infections generated by the virus that is released into the environment during their life cycle.

From (3.4), we can know that the basic reproduction number can be decomposed into two parts, i.e., the direct, , and indirect reproduction number, . Because of the complexity of , defined as the spectral radius of , and hence, in the numerical simulation in Section 4.2.2, is modeled as a constant to be estimated, denoted by , to be a periodic time-varying function constructed by a step function. According to the strong seasonality in both TB incidences and relative infectivity, the periodicity of is considered to be one year, i.e., . We model changes across different months, in other words, we can obtain different values of s in each month. The value of is restricted to be the same within each month.

Threshold dynamics

If , the DFE of model (2.3) is global asymptotically stable.

Consider an auxiliary systemwhich is equivalent towhere and

It follows Lemma 2.2 that there exits a positive -periodic function such that is a solution of (3.5), where . Choose and a small number such that . Then we can get for . By the comparison principle, we havewhere is the transposition of the vector. It follows from that . Since is continuous for all small , we can choose small enough such that . Hence, we get . It follows that as . Hence . By the fourth and fifth equation of model (2.3), we get . This indicates that DFE of model (2.3) is global asymptotically stable.

In the following, we attempt to explore the uniform persistence of model (2.1) when . Define

X and are positively invariant.

From the fifth equation of model (2.3), we can derive thatwhere . Then

By (Smith, 1996, Theorem 4.1.1) as generalized to nonautonomous systems, the irreducibility of the cooperative matriximplies that . Thus, X and are positively invariant. Clearly, is relatively closed in X.

Let be the Poincaré map associated with model (2.3), that is,where is the unique solution of model (2.3) with . Then

From Theorem 2.1, we know that P is a dissipative point on . Thus P admits a global attractor, which attracts every bounded set in . We then introduce the following lemma.

If , there exists a such that, , when , there is where represents distance between and .

Since , Lemma 3.1 implies that . Then we can choose small enough that , where

By the continuity of the solutions with respect to the initial values, for , there exists a such that with , we can obtain

We proceed by contradiction to prove that

If not, we can get

Without loss of generality, we can assume that for all . Then we can get

For any , let , where and , which is the greatest integer less than or equal to . Then,

It follows from (2.5) that there exists that and for .

Consider the following auxiliary systemwhich is equivalent towhere and defined as in (3.7).

It follows from Lemma 2.2 that there exists a positive -periodic function such that is a solution of (3.5), where . Choose and a small number such that , Then we have get for . By the comparison principle, we haveThen as , a contradiction. This completes the proof.

Define

P is uniformly persistent with respect to .

Now we first prove

Noting thatwe only need to prove that

It suffices to prove that for any , we have . If it is not true, there exists an , such that . Thus (3.6) implies thatby replacing the initial time 0 with . Similarly, By (Smith, 1996, Theorem 4.1.1) as generalized to nonautonomous systems, it follows that . where the initial value . Then we have , i.e,

Thus, if , then , which contradicts with . Hence, ,which implies that (3.9) holds. Clearly, is the only fixed point of P and acyclic in . Moreover, Lemma 3.6 implies that is an isolated invariant set in X and , where is the stable set of . By the acyclicity theorem on uniform persistence for maps  citep[Thorem 3.1.1]zhao2003dynamical, it follows that P is uniformly persistent with respect to .

(Zhao et al., 2003, p.18) (Uniformly persistent) Model (2.3) is said to be uniformly persistent if there exists a constant such that any solution with satisfies

If , model (2.3) has at least one positive periodic solution which is uniformly persistent.

It follows from Lemma 3.7 and (Zhao et al., 2003, Thorem 3.1.1) that the solution of model (2.3) is uniformly persistent.

Furthermore, taking advantage of Zhao et al. (2003, Theorem 1.3.6), P has a fixed point . Then, we see that . We further prove that . Suppose not, if , form the last equation of model (2.3), we derive thatwith , where . Therefore, we havewhich yields a contradiction. Hence, and is a positive -periodic solution of model (2.3). This completes the proof.

We next consider the special case of . For simplicity, definewhere .

To find endemic equilibrium, we make the substitution . Then,

Since , By using some algebraic computations, we can obtainwhere x is a positive real root of the following equation:where,and

If , then , and . It follows that Eqn. (3.13) has no positive real root. If , Eq. (3.13) has a unique real root :

Hence model (2.3) has a unique endemic equilibrium with

In a special case of , if , model (2.3) has a unique endemic equilibrium with

TB epidemics in Jiangsu, China via numerical simulations

TB epidemics in Jiangsu, China from 2009 to 2018

The monthly TB incident cases are collected from the Jiangsu provincial center for diseases control and prevention (CDC) (The reported tuberculosis cases in Jiangsu province, 2018). In Fig. 3(a), we show the local TB epidemic from 2009 to 2018. We can observe that a substantially decreasing trend in the TB incidences, dropped from some 4000 cases per month in 2009 to some 2000 cases per month in 2018.

The relative infectivity can be (preliminary) quantified by using the approach in Fine and Clarkson (1982) as well as adopted in Zhao et al. (2018). The relative infectivity can be easily calculated by using the ratio of the number of incidences of time to the number of incidences of time t, i.e., . This Q appears to be a simplified version of quantifying the time-varying (effective) reproduction number by the serial interval approach as studied and implemented in Wallinga and Teunis, 2004, Fraser, 2007, Zhao et al., 2019. The relative infectivity of TB, , is quantified in Fig. 3(b).

In Fig. 3(c), we show the strong seasonality (Fine and Clarkson, 1982) in the TB incidence time series in the annualised epidemic curves. And in Fig. 3(d), we show the annualised relative infectivity to show the seasonality in the TB infectivity across years. We can find that the (relative) infectivity in February appeared to be dominant (or the highest) across different months.

Fitting and estimation results via numerical simulations

Statistical fitting framework

Based on the epidemic model (2.1), we compute the monthly number of reported cases, , of the i-th month (during the study period) aswhere is a constant scaling term for the number of TB cases. In other words, represents a combined effect of the TB symptomatic rate and the reporting rate. Obviously, denotes the theoretical monthly TB cases yielding from model (2.1).

On the other hand, we treat the observed (or reported) number of TB cases, for the i-th month, as a partially observed Markov process (POMP) (King et al., 2016), also know as the hidden Markov model (HMM) from the theoretical number of cases, i.e., in Eqn. (4.1).

We adopt the Poisson-distributed priors for the s such that all s are assumed to follow Poisson distributions according to the theoretical outcomes, i.e., s (Zhao et al., 2018). In other words, the rate of Poisson distribution is a variable depending on , and the observed number of TB cases, , is a random sample from the (predetermined) Poisson distribution in Eqn (4.2). Therefore,

We denote to be the likelihood function of the i-th month, which is the measurement of the “probability” of the observed , given the theoretical number of cases being under the Poisson distribution (Zhao et al., 2018, He et al., 2009, Lin et al., 2018).

Gathering all s, the overall log-likelihood, denoted by l, for the whole TB incidences time series is given in Eqn (4.3).where is the parameter vector to be estimated. The term M denotes the total number of months during the study period, i.e., from 2009 to 2018. We apply the plug-and-play likelihood-based inference framework to estimate the maximum likelihood estimates (MLE) of (He et al., 2009, Lin et al., 2018, Ionides et al., 2006). The profile likelihood approach is implemented to inference the confidence intervals of the model parameters to be estimated (Ionides et al., 2006, Ionides et al., 2017). We use the fixed-time-step Euler-multinomial algorithm (Zhao et al., 2018, Lin et al., 2018) to simulate the epidemic model (2.1).

We consider that there are equivalent birth and death rate (by forcing in model (2.1)), and zero disease-induced mortality rate () in the whole study period, from 2009 to 2018. In this case, the number of the total population, N, is a constant, which, in Jiangsu, slightly changed from 78.1 million in 2009 to 80.5 million in 2018 (The reported tuberculosis cases in Jiangsu province, 2018).

We attempt to find the MLEs of both and within biologically and clinically reasonable ranges by seeking for maximal value(s) of l in (4.3). And this reconstruction approach allows us to project the TB epidemics in an intuitive manner. The projection are conducted by simulating the model (2.1) with MLEs of the parameters up to the end of 2019.

The model simulations are conducted by using the software R (version 3.6.3) (Team, 2013).

Fitting and estimation results

First of all, we show the parameters’ values used for the numerical simulation and sensitivity analysis for model (2.1) in Table 1 .

Table 1

The summary table of model parameters’ values.

parameter	value	unit	statu
γ1-1=γ2-1	1.5	year	fixed
μ^-1	75	per year	fixed
Λ	μN	person	fixed
p	0.05	per year	fixed
q	0.1	per year	fixed
α	1	per case day	assumed
c	0.1	per day	assumed
κ	0.01	unit-free	fixed
ρ	0	per day	assumed
v	1	per day	assumed
β_1(t)=β_2(t)	time-varying	per day	to be estimated
N	8×10^7	person	fixed
S(0)	0.15	unit-free	fixed
T(0)	1×10^-3	unit-free	fixed
T_e(0)	1×10^-3	unit-free	fixed
R(0)	1-S(0)-T(0)-T_e(0)	unit-free	fixed
W(0)	0.04	unit-free	assumed

In Fig. 4 (a)-(b), we show the MLEs of and of each month, respectively. And the annual average is estimated of 1.13, with of 0.35 (see Fig. 5 (a) for details). The is estimated to be strictly less than one. We can find that the estimated annual trends in (Fig. 4(b)) are consistent with the patterns in the preliminary descriptive “infectivity” in Fig. 3(d). We estimated that the in February is larger than one (estimated to be 6.9 in Fig. 5(b)), whereas those in other months are likely to below one.

Fig. 4

The estimation of direct, , and indirect reproduction number, , and the model simulation and projection results. Panel (a) shows the maximal likelihood estimation (MLE) of . Panel (b) shows the MLEs of the of each month. Panel (c) presents the model simulation results, the black dots are the observed number of incidence, the blue line is the model fitting result and the green dashed line is the model projection result. The shading areas represents the 95% credible intervals (CI). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

(a) The profile likelihood of the direct basic reproduction number ; (b) The profile likelihood of the indirect basic reproduction number, in February as an example. The green dots are the random prior samples of different set of parameter values for further simulation purpose. The green curve is the smoothed (by the locally estimated scatterplot smoothing, LOESS) likelihood profile. The horizontal black dashed line is the 95% CI cutoff. The red triangle is the MLE of the parameter of interest. The two vertical red dashed lines indicate the lower and upper bounds of the 95% CI. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In Fig. 4(c), we give the fitting results and projection to 2019. We find that a decreasing trend in the fitting results, which matched the trends in the observed incidence data. The projection is also likely to maintain these trends in 2019.

It is should be noted that the term is considered as a time-varying parameter such that the indirect transmission could also be time-varying. Hence, the reconstructed in Fig. 4(b) is the time-varying reproduction number.

Getting back to the basic reproduction number defined in (3.2), we agree that it should be in an autonomous setting that all model parameters are fixed. In this case, the term in (3.2) is the average of , i.e., annual average. Thus, the basic reproduction number, , in (3.2) is the annual average of the time-varying, which is estimated at 1.13 in Fig. 4.

In addition, in the special case of and , since, in (3.4), and and are all one-to-one mapping, the profile of fitted values of and can thus be directly derived from the estimates in Fig. 5 and other fixed parameters in Table 1. We consider the terms and appear easier to interpret and more biologically meaningful, and thus we choose to show the fitting values of and in Fig. 5.

Since the numerical simulation is conducted with the same complexity in the model structure, we can directly study the goodness-of-fit (in term of the likelihood) and the fitting errors. The goodness-of-fit and the error term analyses are demonstrated in Fig. 6 . We can find that the fitted values of the TB incidences are in line with the observations with the error terms (largely) following in a Normally distribution (Fig. 6(a)-(b)). And the mean percentage error is close to zero (Fig. 6(c)).

Fig. 6

The matching between the observed and fitted values and the distribution of the fitting error terms, e.g., the differences between observed and fitted values. Panel (a) shows the (normalised) observations against the fits (dots), and the diagonal line represents the “” line. Panel (b) shows the Normal quantile–quantile (QQ) plot of the fitting errors. Panel (c) shows the distribution of the percentage errors, i.e., the ratios of the errors over the observations. Panel (d) shows the distribution of the errors.

Sensitivity analysis and the trend change in the TB epidemic

Following (Zhao et al., 2018, Zhao et al., 2018, Musa et al., 2019, Gao et al., 2016, Tang et al., 2016, Tang et al., 2016), we adopt the partial ranked correlation coefficient (PRCC) for the sensitivity analysis between the model outcomes and the parameters. The PRCCs of the , infection attack rate (IAR) and the environmental contamination level of the model (2.1) are estimated. The sensitivity analysis results are in Fig. 7 , and suggest that most of the model parameters are significantly associated with the TB infectivity, IAR and the environmental contamination, which should be given priorities in controlling the TB epidemics.

Fig. 7

The partial rank correlation coefficients (PRCC) of basic reproduction number in panel (a), the infection attack rate (IAR) in panel (b) and the level of the environmental contamination in panel (c) against the model parameters. The denotes the initial susceptible ratio. The denotes the initial environmental contamination level. The dots are the estimated PRCCs, and the bars represent the 95% CIs. The ranges of model parameters are based on the values in Table 1 having a random perturbation with a coefficient of variation of 0.2.

Based on the results of the sensitivity analysis above, we conduct the numerical simulations to present the changing dynamics of the TB epidemics and the environmental contamination levels with changes in the epidemiological parameters. Fig. 8 (a) and (b) show the trend changes in the TB epidemic with changes in the and , respectively. Fig. 8(c) and (d) show the trend changes in the environmental contamination levels with changes in the parameters and c, respectively.

Fig. 8

The numerical simulation results with the changes in epidemiological parameters. The panels (a) and (b) are the numbers of TB cases. The panels (c) and (d) are the levels of the environmental contamination. In panels (a) and (b), the blue lines are the same main results, by using the MLEs, as in Fig. 4(c). In panels (c) and (d), the black lines are the same main results as in Fig. 4(c). Except for those indicated in the figure legends, all other model parameters and initial conditions are the same as in Table 1.

Concluding remarks

The main focus of this study is to investigate the effects of the contaminated environments on the TB transmission dynamics in Jiangsu, China analytically and numerically. Mathematically, we define the basic reproduction number (cf. (3.2  ), and prove that can be used to govern the threshold dynamics of the model: if , the unique DFE is globally asymptotic stable (cf. Theorem 3.4); while , there is at least one positive periodic solution and TB will persist uniformly (Theorem 3.9). Epidemiologically, we show that the cost of the contaminated environments affect the transmission dynamics of TB in Jiangsu, China in the following aspects:

(i)  Based on the monthly TB incident cases counted by the Jiangsu CDC (cf. Fig. 3(a) ), the TB incidence time series has strong seasonality (cf. Fig. 3(c) ). And the annual average is , then from Theorem 3.9, we can conclude that the TB in Jiangsu persists under current circumstances. That is, the TB becomes an endemic disease and will persist in Jiangsu for a long time. And there is a long way to go to achieve the world-wide goal towards elimination of tuberculosis by 2050.

(ii)  The annualised relative infectivity (cf. Fig. 3(d) ) shows that the relative infectivity in February is dominant across different months. This is consistent with the estimation of that in February is 6.9 (cf. Fig. 5(b) ), whereas those in other months are likely to below one. This phenomenon seems to be the first reported case. For one of the possible explanations, we conjecture that this is induced by the travel or migration with the winter vacation and the Spring Festival, commonly started since mid-February of each year and lasted for about 20 days. During this period, a large number of people working or living outside of Jiangsu will return to their hometown, and the local CDC will increase the screening of tuberculosis for returning people. Moreover, due to population movements and increased exposure rates, the risk of tuberculosis transmission has increased significantly, and the resulting lagging effect will increase the number of cases in the next month (see Fig. 3(c) ). Although this has not yet been formally verified in public health filed, it is desirable in future studies. From the model fitting side, both of the surveillance reporting and the TB infectivity have similar effects on the number of TB incidences. This means it is difficult to disentangle the solo effect of both factor as the same time based on our model framework. We choose to set the surveillance reporting efforts, in term of in Eqn (4.1), as a constant during the entire study period, i.e., 2009–2018, and estimate the TB infectivity, in term of , as a time-varying function. If both and were set to be time-varying in the fitting procedure, the potential over-fitting problem as well as the estimation biases in both and would probably occur. This is largely due to the effects of and cannot be disentangled and thus not independent in our model structure. Although lack of supporting information or data, our model is still capable to capture the long-term TB epidemic in Jiangsu. We remark that more detailed information on quantifying the local TB surveillance and the changing dynamics of the de jure population and floating population would be very helpful to identifying more accurate TB infectivity estimates.

(iii)  From the numerical results in Fig. 5, we find that merely controlling the changing dynamics of the TB indirect transmission, in term of the , appears sufficient to successfully capture the long-term patterns in TB epidemics in Jiangsu. With estimated strictly less than one, we remark that the TB epidemics are likely to be controlled, in term of , by effectively controlling the indirect transmission path. Also, Fig. 8 indicates that the control measures reduce the or could effectively decrease the number of TB cases. This could be achieved by providing timely and effective treatment and maintaining a low environment contamination level. We further find that a lower virus shedding rate, , and a higher environmental clearance rate, c, lead to low level of environment contamination (Fig. 8) and (Fig. 7), which could control the TB epidemic efficiently.

(iv) From Fig. 7, the TB transmissibility and number of cases are positively associated with the effective transmission rates and , as well as the virus shedding rate . The effective control efforts are suggested to focus on reducing the and . We also find that the recovery rates, and , and environmental clearance rate, c, are negatively associated with the TB transmissibility and number of cases. Thus, increasing and are also likely to control the TB epidemics.

It is worthy to note that, increasing number of evidences support that the respiratory infections are primarily transmitted between people through respiratory droplets and contact routes World Health, 2020, Liu et al., 2020, Li et al., 2020, Huang et al., 2020). Recently, Gao et al. (2020) investigated the relative contributions of different transmission routes to a multi-route transmitted respiratory infection and found that all transmission routes can dominate the total transmission risk under different scenarios. In the present paper, we model the transmissions of the TB into two categories, one is direct, and the other is indirect, which is measure by representing the effects of the contaminated environments. Our model (2.1) provides a straightforward method to evaluate the transmission efficiency of different transmission routes of TB.

In addition, the progress in controlling TB is, however, currently influenced by some major factors, such as multidrug resistant (MDR) (Dodd et al., 2016, Knight et al., 2019, Liu et al., 2019), ambient particulate air pollution (Liu et al., 2019, Peng et al., 2017, Liu et al., 2019), etc. The effects of MDR or air pollution on the transmission dynamics of TB in Jiangsu, China will be desirable in our future studies.

CRediT authorship contribution statement

Yongli Cai: Modelling, mathematical analysis, writing. Shi Zhao: Modelling, numerical analysis, writing. Yun Niu: Modelling, numerical analysis. Zhihang Peng: Modelling, mathematical analysis. Kai Wang: Modelling, mathematical analysis. Daihai He: Modelling, numerical analysis. Weiming Wang: Modelling, mathematical analysis, numerical analysis, Writing-Reviewing and Editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.