Impact of variability of reproductive ageing and rate on childhood infectious disease prevention and control: insights from stage-structured population models

1.   Introduction

Family planning involves consideration of the number (including the choice of zero) and spacing of children a family wishes to have. A number of factors can impact the family planning at the individual family level. At the population level the family planning and reproduction strategy including fertility, birth age and spacing of children, may be heavily influenced by economical conditions and societal resources which can be weighted heavily by the age-distribution of the entire population. In developing countries, policies like subsidizing education raise the earning power of women and the opportunity cost of having children, consequently lowers fertility ^[1]. Access to contraceptives may also yield lower fertility rates. In developed countries, the proportion of retired people is increasing, adding burden on the workforce population to support pensions and social programs. Increasing high skill migration may be an effective way to increase the return to education leading to lower fertility and a greater supply of highly skilled individuals ^[1], thus address the aging population problem.

A well-known example of family planning and age-distribution of population being significantly regulated by political and social-economic consideration is the China's one-child policy implemented for many years. In 1973, the Chinese government issued voluntary guidelines on fertility control to encourage later marriage, longer spacing between births, and fewer births overall ^[2,3]. In 1981, China's National Family Planning Commission proposed a population control policy advocating one child per couple, which was moderated in 1984, allowing most rural families a second child ^[2,3,4]. In 2002, the policy was incorporated in the Population and Family Planning Law, at the same time, a second child was permitted in some provinces if both husband and wife were from single-child families ^[2,3,5]. In 2013, the policy was relaxed to allow a second child if either spouse was from a single-child family ^[3,6]. In October of 2015, the Chinese government announced a two-child policy, effective from January 1 of 2016 ^[7]. The new policy that allows each couple have two children was proposed in order to help address the population aging issue. It was reported that, starting from May 2018, Chinese authorities were in the process of ending the population control policies ^[8].

A consequence of this recent change of the centralized population control policy after a long-term implementation of one-child per family policy is the obvious increasing of the family size, and substantial heterogeneity of the reproduction age and the spacing between the first and second child in those families with two children. This generates new close contact patterns in household and community level and thus any issue relevant to these contact patterns must be revisited. The control and prevention of childhood infectious diseases preventable by vaccine, such as pertussis, is one of these critical public health issues. Taking pertussis as an example, this childhood disease can be fatal in infants but infection can be prevented in other age groups with an effective vaccine. Pertussis vaccines wane over time, so those children who are expected to have younger siblings need to take a booster vaccine if (I) the prevalence of disease in the older age group and/or groups (recalling the potential heterogeneity of spacing between two children since females in multiple age groups may consider to give birth) is expected to be high; and (II) vaccine waning make this group of these groups less protected and more susceptible to the disease.

To the best of my knowledge, there is no study on impact of family planning and the scale of density-regulated birth rate on the long-term population demographic distribution and childhood disease dynamics. However, there are a few studies which imply the impact of demographic change on infectious disease dynamics. These work include studies on demographic transition and the dynamics of measles in China ^[9], the influence of demographic change on spread of infectious diseases ^[10], the impact of demographic transition on rubella transmission dynamics in China ^[11], the effects of demographic change and immigration on infectious diseases in Italy ^[12], the effects of demographic change on disease transmission and vaccine impact in a household structured population ^[13], the dynamical consequences of demographic change in a model of disease transmission ^[14] and the impact of demographic change on the estimated future burden of hepatitis B and seasonal influenza in the Netherlands ^[15].

This series of studies is dedicated to developing mathematical frameworks and analyses to examine the patterns of childhood infectious disease transmission, to identify prevalence of disease in different age groups, when female in multiple age intervals are giving birth to the second children. In this first paper of the series, we start with a simple stage-structured disease transmission model, and study the impact of family planning and the scale of density-regulated birth rate on the long-term population demographic distribution and infant disease incidence prevalence. We conduct our analyses by varying three parameters: birth rate, reproduction age interval(s), and the scale of the sub-population density regulation.

In particular, we introduce a multi-stage ($m$-stages) stratified model, where the population is divided by age into $m$ groups (stages) with the $i$-th age stage spanning the age interval of length $\tau_{i}$. Assume that females in the age groups, $k, k+1, \cdots, l$-th groups give birth. For each group, we have the classical SIS epidemic-model, where the population is divided into the susceptible and the infectious. We consider the situation where the infectious period is much shorter than the period of each age stage. Therefore, the infectious individuals in the $i$-th age group $I_{i}$ go back to the susceptible class $S_{i}$ before advancing to the $(i+1)$-th age group. The flowcharts of the demographic model and epidemiological model are shown in Figures 1 and 2 respectively.

Let $S_{i}(t)$ be the population of the susceptible of the $i$th age group, $I_{i}(t)$ be the population of the infectious of the $i$th age group. $N_{i}(t)$ denotes the total population of the $i$th age group at time $t$. The death rate of the $i$th age group is given by $\mu_{i}$; $\sigma$ is the recover rate; the birth rate of the $i$th productive group is a nonlinear function $b_{i}(N_{i})$; $\beta_{ij}$ is the transmission rate of the disease from stage $j$ to stage $i$. The age-stratified epidemiological model is given by the following equations:

The demographic model is given by

2.   The demographic model

In this section, we investigate the dynamics of the demographic model by studying the stability of equilibrium.

Linearizing (1.2) at the zero equilibrium gives

Let $\lambda$ be eigenvalue of the linear system (2.1). By calculation, the characteristic equation is given by

Now we make the following assumption:

(A1) The birth function takes the form $b_{i}(x) = p_{i}xq(x)$ where $q(x)$ is a non-negative monotone decreasing function.

Note that $p_{i}$ is the maximal number of children a female in age group $i$ could give per unit time, $q(x)$ is the function which implies the restriction of resources, so assumption (A1) reflects the ecological consideration that the reproduction is linear in $x$ only for small densities and decreases as a consequence of intra specific competition. For example, one well known birth function which takes the form in assumption (A1) is the Ricker function $b(x) = pxe^{-qx}$.

With this assumption, we have the following theorem on local stability of the zero equilibrium.

Theorem 1. Under assumptions (A1), if $\sum_{i = k}^{l}p_{i}q(0)\frac{e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1-e^{-\mu_{i}\tau_{i\;}})}{\mu_{i}} > 1$, then the zero equilibrium is unstable; if $\sum_{i = k}^{l}p_{i}q(0)\frac{e^{-\sum_{j = 1\;}^{i-1\;}\mu_{j}\tau_{j}}(1-e^{-\mu_{i}\tau_{i}}\;)}{\mu_{i}} < 1$, the zero equilibrium is stable.

Proof. Let $G(\lambda): = \sum_{i = k}^{l}b_{i}'(0)g_{i}(\lambda)$, where $g_{i}(\lambda) = \frac{(1-e^{-(\lambda+\mu_{i})\tau_{i}})e^{-\lambda\sum_{j = 1\;}^{i-1}\tau_{j}-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}}{\mu_{i}+\lambda}$. Then the characteristic equation (2.2) can be written as $G(\lambda) = 1$. Calculating the derivative gives $g_{i}'(\lambda) < 0$ on $(-\mu_{i}, \infty)$. Furthermore, $g_{i}(\lambda)\rightarrow 0$ as $\lambda\rightarrow\infty$ and $g_{i}(\lambda)\rightarrow +\infty$ as $\lambda\rightarrow -\mu_{i}$. From assumption (A1), we have $b_{i}'(0) = p_{i}q(0) > 0$. Therefore, $G(\lambda)$ is monotone decreasing on $(-\mu, \infty)$ where $\mu = \min\{\mu_{i}, i = k, k+1, ..., l\}$. Moreover, $G(\lambda)\rightarrow 0$ as $\lambda\rightarrow\infty$ and $G(\lambda)\rightarrow +\infty$ as $\lambda\rightarrow -\mu$.

If $G(0) = \sum_{i = k}^{l}p_{i}q(0)\frac{e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1-e^{-\mu_{i}\tau_{i}})}{\mu_{i}} > 1$, $G(\lambda) > 1$, since $G(\lambda)$ is monotone decreasing on $(-\mu, \infty)$ and $\lim_{\lambda\rightarrow\infty}G(\lambda) = 0$, the characteristic equation has a positive real root. So the zero equilibrium is unstable. If $G(0) = \sum_{i = k}^{l}p_{i}q(0)\frac{e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1-e^{-\mu_{i}\tau_{i}})}{\mu_{i}} < 1$, $G(\lambda) < 1$, since $G(\lambda)$ is monotone decreasing on $(-\mu, \infty)$, $\lim_{\lambda\rightarrow\infty}G(\lambda) = 0$ and $\lim_{\lambda\rightarrow -\mu}G(\lambda) = +\infty$, the characteristic equation has a positive real root. So the zero equilibrium is stable.

Suppose that there is a positive equilibrium $(N_{1}^{*}, N_{2}^{*}, ..., N_{m}^{*})$, then we have

From Eq (2.3) we derive

The conditions for existence of this positive equilibrium is given in the following theorem.

Theorem 2. Under assumption (A1) with $\lim_{x\rightarrow\infty}q(x) = 0$, the positive equilibrium $(N_{1}^{*}, N_{2}^{*}, ..., N_{m}^{*})$ exists and is unique if $\sum_{i = k}^{l}p_{i}q(0)\frac{e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1-e^{-\mu_{i}\tau_{i}})}{\mu_{i}} > 1$.

Proof. From Eq (2.4),

for $k\leq h\leq l$. Equation (2.5) and the first equation in (2.4) imply that

From assumption (A1), $b_{i}(x) = p_{i}xq(x)$, so Eq (2.6) becomes

So the positive equilibrium exists if there exists a positive $N_{1}^{*}$ such that Eq (2.7) holds. Now let $G(x) = \sum_{i = k}^{l}p_{i}\frac{e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1-e^{-\mu_{i}\tau_{i}})}{\mu_{i}}q(\frac{\mu_{1}e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1-e^{-\mu_{i}\tau_{i}})}{\mu_{i}(1-e^{-\mu_{1}\tau_{1}})}x)$. Since $q(x)$ is monotone decreasing with respect to $x$, $G(x)$ is monotone decreasing function. Furthermore, $\lim_{x\rightarrow\infty}G(x) = 0$. So $G(x) = 1$ has a unique positive solution if and only if $G(0) > 1$, i.e., $\sum_{i = k}^{l}p_{i}q(0)\frac{e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1-e^{-\mu_{i}\tau_{i}})}{\mu_{i}} > 1$.

Note that Theorems 1 and 2 imply that the positive equilibrium exists and is unique if and only if the zero equilibrium is unstable.

For the next, we study stability of this positive equilibrium.

We denote by $C_{m}^{+}$ the non-negative cone of the Banach space of continuous functions $C_{m} = \{\varphi = (\varphi_{1}, \varphi_{2}, ..., \varphi_{m}):[-r, 0]\rightarrow \mathbb{R}^{m} continuous\}$, where $r = \max\{\tau_{1}, \tau_{2}, ..., \tau_{m}\}$, i.e. $C_{m}^{+} = \{\varphi\in C_{m}: \varphi_{i}(\theta)\geq 0 \quad for\quad \theta\in[-r, 0], i = 0, 1, 2, ..., m\}$. By using the method of steps, it can be shown that for each $\varphi\in C_{m}^{+}$, there is a unique solution of (1.2) $\pi(\varphi, t) = (N_{1}(\varphi, t), N_{2}(\varphi, t), ..., N_{m}(\varphi, t))\in\mathbb{R}_{m}^{+}$ through $\varphi$ that is well defined and satisfies $\pi(\varphi; .)|_{[-r, 0]} = \varphi$.

In fact, by taking integral and making substitutions, system (1.2) can be written as

In what follows, we give a preliminary result, then we give a theorem on global stability of the positive equilibrium.

Lemma 1. Under assumption (A1), if the birth functions $b_{i}(x)$ are bounded for $i = k, k+1, ..., l$, for every $\varphi\in C_{m}^{+}$ with $\varphi_{i}(0) > 0$, $i = 1, 2, ..., m$, the solution $\pi(\varphi; t)$ of (1.2) is bounded above for $t > 0$.

Proof. Let $N(t) = \sum_{i = 1}^{m}N_{i}(t)$. By adding up the $m$ equations in Eq (1.2), we obtain

where $\mu$ is the smallest death rate in the $m$ age groups, i.e., $\mu = \min\{\mu_{i}, i = 1, 2, ..., m\}$. Since the birth functions $b_{i}(x)$ are bounded for $i = k, k+1, ..., l$, there are $M_{i}$ for $i = 1, 2, ..., m$ such that $b_{i}(N_{i}(t))\leq M_{i}$. Let $M = \sum_{i = k}^{l}M_{i}$, then $\frac{dN}{dt}\leq M-\mu N(t)$, which means that $\frac{dN}{dt} < 0$ when $N > \frac{M}{\mu}$. So $N$ is bounded, i.e., there is $\bar N$ such that $N(t)\leq \bar N$ for $t\geq 0$. Therefore, $N_{i}(t)\leq\bar N$ for $t\geq 0$ for $i = 1, 2, ..., m$. The solution $\pi(\varphi; t)$ is bounded for $t > 0$.

Theorem 3. Under assumption (A1) with $\lim_{x\rightarrow\infty}q(x) = 0$, if $\sum_{i = k}^{l}p_{i}q(0)\frac{e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1-e^{-\mu_{i}\tau_{i}})}{\mu_{i}} > 1$ and $\sum_{i = k}^{l}\frac{\left\vert b_{i}'(N_{i}^{*})\right\vert e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1+e^{-\mu_{i}\tau_{i}})}{\mu_{i}} < 1$, then the positive equilibrium $(N_{1}^{*}, N_{2}^{*}, ..., N_{m}^{*})$ is locally stable.

Proof. The linearized equations at the endemic equilirbium $(N_{1}^{*}, N_{2}^{*}, ..., N_{m}^{*})$ of system (1.2) is given by

Let $\lambda$ be eigenvalue of the linear system (2.9). By calculation, the characteristic equation is given by

Suppose that the characteristic equation (2.10) has an eigenvalue with non-negative real part, i.e., there exits $\lambda = x+iy$ such that $x\geq0$, then

Therefore, Eq (2.10) and inequality (2.11) indicate that

which contradicts with the assumption that $\sum_{i = k}^{l}\frac{\left\vert b_{i}'(N_{i}^{*})\right\vert e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1+e^{-\mu_{i}\tau_{i}})}{\mu_{i}} < 1$. So the characteristic equation (2.10) has no eigenvalue with non-negative real part, the positive equilibrium $(N_{1}^{*}, N_{2}^{*}, ..., N_{m}^{*})$ is locally stable.

Theorem 4. Under assumption (A1) with $\lim_{x\rightarrow\infty}q(x) = 0$ and the birth functions $b_{i}(x)$ bounded for $i = k, k+1, ..., l$, assume that $\mu_{h} = \mu$ and $\tau_{h} = \tau$ for some $\mu > 0$, $\tau > 0$ and all $k\leq h\leq l$. If $\sum_{i = k}^{l}\frac{\left\vert b_{i}'(N_{i}^{*})\right\vert e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1+e^{-\mu_{i}\tau_{i}})}{\mu_{i}} < 1$, then the positive equilibrium $(N_{1}^{*}, N_{2}^{*}, ..., N_{m}^{*})$ is globally stable. i.e., $\lim_{t\rightarrow\infty}\pi(\varphi; t) = (N_{1}^{*}, N_{2}^{*}, ..., N_{m}^{*})$ for $\varphi\in C_{m}^{+}$ with $\varphi_{i}(0) > 0$.

Proof. Let $\{N_{i}(t)\}$ be a solution of Eq (1.2). Since it's bounded, we can define

Let $h$ be such that $k\leq h\leq l$, i.e. $N_{h}$ is a productive group. There exists a sequence $\{t_{n}\}$ and a sequence $\{s_{n}\}$ such that $\lim_{n\rightarrow\infty}N_{h}(t_{n}) = \gamma_{h}$ and $\lim_{n\rightarrow\infty}N_{h}(s_{n}) = \delta_{h}$. So there exists some $\epsilon > 0$ such that $\delta_{i}-\epsilon < N_{i}(t_{n}) < \gamma_{i}+\epsilon$ and $\delta_{i}-\epsilon < N_{i}(s_{n}) < \gamma_{i}+\epsilon$ for $n$ large enough for all $k\leq i\leq l$.

From the integrated equation (2.8),

Let $n\rightarrow\infty$ and $\epsilon\rightarrow 0$, the inequality becomes

Now let $A: = \sum_{i = k}^{l}p_{i}\gamma_{i}q(\delta_{i})$, Eq (2.13) implies that

which further implies that

Let $B: = \sum_{i = k}^{l}p_{i}N_{i}^{*}q(N_{i}^{*})$, from Eq (2.4)

so

From assumption (A1), $q(x)$ is monotone decreasing, then Eqs (2.14) and (2.15) imply that there exists $h_{1}$ such that $k\leq h_{1}\leq l$ and $\delta_{h_{1}}\leq N_{h_{1}}^{*}$. From the integrated equation (2.8),

Following similar calculation as inequalities (2.13) and (2.14), we obtain

From assumption (A1), $q(x)$ is monotone decreasing, then Eqs (2.15) and (2.16) imply that there exists $h_{2}$ such that $k\leq h_{2}\leq l$ and $\gamma_{h_{2}}\geq N_{h_{2}}^{*}$.

By substituting variables, from Eqs (2.4) and (2.8) we obtain

and

Since $\mu_{h} = \mu$ and $\tau_{h} = \tau$ for $k\leq h\leq l$, Eqs (2.17) and (2.18) imply that

for $k\leq h\leq l$.

Scenario 1: $\gamma_{k} = N_{k}^{*}$

In this scenario, since the positive equilibrium is locally stable, $\delta_{k} = \gamma_{k} = N_{k}^{*}$, and $\lim_{t\rightarrow\infty}N_{k}(t) = N_{k}^{*}$. From Eq (2.19), $\lim_{t\rightarrow\infty}N_{h}(t) = N_{h}^{*}$ for all $k\leq h\leq l$.

Scenario 2: $\gamma_{k} > N_{k}^{*}$

In this scenario, since the positive equilibrium is locally stable, $\delta_{k} > N_{k}^{*}$, then $N_{k}(t) > N_{k}^{*}$ for $t$ large enough. From Eq (2.19), $N_{h_{1}}(t) > N_{h_{1}}^{*}$ for $t$ large enough, which means that $\delta_{h_{1}}\geq N_{h_{1}}^{*}$. Since we have $\delta_{h_{1}}\leq N_{h_{1}}^{*}$ from previous discussion, $\delta_{h_{1}} = N_{h_{1}}^{*}$. Therefore, $\delta_{h_{1}} = \gamma_{h_{1}} = N_{h_{1}}^{*}$, i.e., $\lim_{t\rightarrow\infty}N_{h_{1}}(t) = N_{h_{1}}^{*}$. From Eq (2.19), $\lim_{t\rightarrow\infty}N_{h}(t) = N_{h}^{*}$ for all $k\leq h\leq l$.

Scenario 3: $\gamma_{k} < N_{k}^{*}$

Then $\delta_{k} < N_{k}^{*}$ and $N_{k}(t) < N_{k}^{*}$ for $t$ large enough. From Eq (2.19), $N_{h_{2}}(t) < N_{h_{2}}^{*}$ for $t$ large enough, which implies that $\gamma_{h_{2}}\leq N_{h_{2}}^{*}$. Since we have $\gamma_{h_{2}}\geq N_{h_{2}}^{*}$ from previous discussion, $\gamma_{h_{2}} = N_{h_{2}}^{*}$. i.e., $\lim_{t\rightarrow\infty}N_{h_{2}}(t) = N_{h_{2}}^{*}$. From Eq (2.19), $\lim_{t\rightarrow\infty}N_{h}(t) = N_{h}^{*}$ for all $k\leq h\leq l$.

From discussion above, $\lim_{t\rightarrow\infty}N_{h}(t) = N_{h}^{*}$ for all $k\leq h\leq l$. Then by the integral equation (2.8), $\lim_{t\rightarrow\infty}N_{h}(t) = N_{h}^{*}$ for all $1\leq h\leq m$.

3.   The epidemic model

Now we focus on the epidemic model, which is an ODE system given by

Note that Eq (3.1) is derived from the last equation in (1.1) by replacing $S_{h}(t)$ by $N_{h}(t)-I_{h}(t)$.

Suppose that the population has reached the positive equilibrium, then system (3.1) is given by

Let $\mathbf{I} = (I_{1}, I_{2}, ..., I_{m})^{T}$, the flow of the solution of system (3.2) with initial value $\mathbf{I^{0}} = (I_{1}^{0}, I_{2}^{0}, ..., I_{m}^{0})$ is given by $\phi_{t}(\mathbf{I^{0}})$. Let $V = (0, N_{1}^{*})\times(0, N_{2}^{*})\times...\times(0, N_{m}^{*})$. We have the following conclusion

Theorem 5. If $\mathbf{I}_{i}^{0}\in(0, N_{i}^{*})$, then $\phi_{ti}(\mathbf{I^{0}})\in(0, N_{i}^{*})$, i.e., $V$ is invariant under the flow $\phi_{t}$.

Proof. Suppose that there is a smallest $t_{0}$ such that there is $j\in\{1, 2, ..., m\}$ such that $I_{j}(t_{0}) = 0$. Since $\mathbf{I}(t)$ is not constant 0, we have $\frac{dI_{j}(t)}{dt}|_{t = t_{0}} = \sum_{i = 1}^{m}\beta_{hi}(N_{h}^{*}-I_{h}(t_{0}))I_{i}(t_{0}) > 0$. On the other hand, since $t_{0}$ is the smallest s.t. $I_{j}(t_{0}) = 0$, $\frac{dI_{j}(t)}{dt}|_{t = t_{0}} = \lim_{\epsilon\rightarrow0}\frac{I_{j}(t_{0}-\epsilon)-I_{j}(t_{0})}{-\epsilon}\leq0$, which is a contradiction. So $I_{j}(t) > 0$ for all $1\leq j\leq m$ and all $t > 0$ with initial value in $V$.

Similarly, suppose that there is a smallest $t_{0}$ such that there is $j\in\{1, 2, ..., m\}$ such that $I_{j}(t_{0}) = N_{j}^{*}$, then from Eq (3.2), $\frac{dI_{j}}{dt}|_{t = t_{0}} = -\sigma I_{h}(t_{0})-\mu_{h}I_{h}(t_{0}) < 0$. On the other hand, since $t_{0}$ is the smallest s.t.$I_{j}(t_{0}) = N_{j}^{*}$, $\frac{dI_{j}(t)}{dt}|_{t = t_{0}} = \lim_{\epsilon\rightarrow0}\frac{I_{j}(t_{0}-\epsilon)-I_{j}(t_{0})}{-\epsilon}\geq0$, which is a contradiction. So $I_{j}(t) < 0$ for all $1\leq j\leq m$ and all $t > 0$ with initial value in $V$.

3.1.   Disease-free equilibrium of the epidemic model

Equation (3.2) has $(0, 0, ..., 0)$ as the disease-free equilibrium. Linearization around this equilibrium gives the following linear system

for $1\leq h\leq m$ Following the method in ^[16], we get

and

Then we have

In particular, if we assume that $\beta_{ij} = \alpha_{i}\lambda_{j}$, it can be calculated from induction that the characteristic equation of $FV^{-1}$ is given by $\lambda^{m-1}(\lambda-\sum_{i = 1}^{m}\frac{\alpha_{i}\lambda_{i}N_{i}^{*}}{\sigma+\mu_{i}}) = 0$, then

The following theorem follows

Theorem 6. Assume that $\beta_{ij} = \alpha_{i}\lambda_{j}$. Then when $\sum_{i = 1}^{m}\frac{\alpha_{i}\lambda_{i}N_{i}^{*}}{\sigma+\mu_{i}} < 1$, the disease-free equilibrium of system (3.2) is stable; when $\sum_{i = 1}^{m}\frac{\alpha_{i}\lambda_{i}N_{i}^{*}}{\sigma+\mu_{i}} > 1$, the disease-free equilibrium of system (3.2) is unstable.

For the assumption $\beta_{ij} = \alpha_{i}\lambda_{j}$, if we assume that the population is homogeneously mixed, $\alpha_{i}$ can be interpreted as susceptibility of age group $i$ and $\lambda_{j}$ can be interpreted as infectivity of age group $j$.

3.2.   Endemic equilibrium of the epidemic model

Suppose system (3.2) has a nontrivial equilibrium $(I_{1}^{*}, I_{2}^{*}, ..., I_{m}^{*})$, by plugging in the Eq (3.2) we derive

for $1\leq h\leq m$. In particular, if we assume that $\beta_{ij} = \alpha_{i}\lambda_{j}$, then Eq (3.4) can be written as

for $1\leq h\leq m$. From Eq (3.5) we get

for $1\leq h\leq m$. and

Now let $M_{i} = \frac{N_{i}^{*}}{I_{i}^{*}}-1$ and $l_{i} = \frac{\sigma+\mu_{i}}{\alpha_{i}}$, then

Plugging $I_{i}^{*} = \frac{N_{i}^{*}}{M_{i}+1}$ into Eq (3.6) with $h = 1$, we have

It follows that

Now define $G:(0, N_{1}^{*})\longrightarrow\mathbb{R}$ by

where $M = \frac{N_{1}^{*}}{x}-1$. Then $I_{1}^{*}$ is a solution of $G(x) = 1$.

Note that $G(x)$ is monotone non-increasing with respect to $x$, $G(x)\longrightarrow 0$ as $x\longrightarrow N_{1}^{*}$ and $G(x)\longrightarrow \sum_{i = 1}^{m}\frac{\lambda_{i}N_{i}^{*}}{l_{i}}$ as $x\longrightarrow 0$. So $G(x) = 1$ has a solution in $(0, N_{1}^{*})$ if and only if $\sum_{i = 1}^{m}\frac{\lambda_{i}N_{i}^{*}}{l_{i}} > 1$, and the solution is unique by monotonicity of $G(x)$.

Note that $R_{0} = \sum_{i = 1}^{m}\frac{\lambda_{i}N_{i}^{*}}{l_{i}}$. We conclude that

Theorem 7. Assume that $\beta_{ij} = \alpha_{i}\lambda_{j}$. The endemic equilibrium of system (3.2) exists and is unique if and only if $R_{0} = \sum_{i = 1}^{m}\frac{\lambda_{i}N_{i}^{*}}{l_{i}} > 1$.

Now we state the following theorem on global stability of the endemic equilibrium.

Theorem 8. Assume that $\beta_{ij} = \alpha_{i}\lambda_{j}$.Let $V = (0, N_{1}^{*})\times(0, N_{2}^{*})\times...\times(0, N_{m}^{*})$. If $R_{0} > 1$, then the endemic equilibrium $\{I_{i}^{*}\}$ of system (3.2) attracts all the forward orbits going through $V$.

Proof. Let $F:V\longrightarrow\mathbb{R}^{m}$ be defined by $F_{h}(I_{1}, I_{2}, ..., I_{m}) = \sum_{i = 1}^{m}\beta_{hi}(N_{h}^{*}-I_{h})I_{i}-\sigma I_{h}-\mu_{h}I_{h}$.

It suffices to prove that

(H1) System (3.2) is cooperative, i.e., $\frac{\partial F_{h}}{\partial I_{j}}\geq0$ for $h\neq j$.

(H2) $F$ is irreducible in the sense that the matrix $[\frac{\partial F_{h}}{\partial I_{j}}]$ is irreducible.

(H3) Solutions of Eq (3.2) with initial value $(I_{1}^{0}, I_{2}^{0}, ..., I_{m}^{0})$ such that $\left\vert I_{h}^{0}\right\vert\leq N_{h}^{*}$ are bounded.

Then by Theorems 1.5 and 2.4 in ^[17], (H1) and (H2) imply that system (3.2) doesn't have a non-constant periodic solution. By Theorem 1.1 in ^[17], (H1) and (H2) also imply that the solution flows of Eq (3.2) going through $V$ have positive derivatives. Then by Theorem 4.1 in ^[17] and (H3) we conclude that almost all forward orbits of $V$ converge to the endemic equilibrium $\{I_{i}^{*}\}$.

(H1) $\frac{\partial F_{h}}{\partial I_{j}} = \beta_{hj}(N_{h}^{*}-I_{h}) > 0$ for $h\neq j$.

(H2) $\frac{\partial F_{h}}{\partial I_{h}} = \beta_{hh}(N_{h}^{*}-I_{h})-\sum_{i = 1}^{m}\beta_{hi}I_{i}-\sigma-\mu_{h}$. Now let $A = [\frac{\partial F_{h}}{\partial I_{j}}]$. $A_{hj} = \frac{\partial F_{h}}{\partial I_{j}} > 0$ for $h\neq j$ by (H1). Suppose there is $1\leq h\leq m$ such that $A_{hh}\leq0$, $A_{hh}^{2} = \sum_{i = 1}^{m}A_{hi}A_{ih} > 0$. Therefore, for each pair of indices $h$ and $j$, there exists a natural number $n$ such that $A_{hj}^{n}$ is positive, which implies that the matrix $A$ is irreducible.

(H3) It can be derived directly from Theorem 5.

Now if we look back on the original epidemic model (3.1), we have the following Theorem from Theorem 8.

Theorem 9. Under assumption (A1) with $\lim_{x\rightarrow\infty}q(x) = 0$, assume that $\beta_{ij} = \alpha_{i}\lambda_{j}$ and $\mu_{h} = \mu$, $\tau_{h} = \tau$ for some $\mu > 0$, $\tau > 0$ and all $k\leq h\leq l$.Let $V = (0, N_{1}^{*})\times(0, N_{2}^{*})\times...\times(0, N_{m}^{*})$. If $\sum_{i = k}^{l}\frac{\left\vert b_{i}'(N_{i}^{*})\right\vert e^{-\sum_{j = 1\;}^{i-1}\mu_{j}\tau_{j}}(1+e^{-\mu_{i}\tau_{i}})}{\mu_{i}} < 1$ and $R_{0} > 1$, then $\{I_{i}^{*}\}$ attracts all the forward orbits going through $V$ in system (3.1).

Proof. From Theorem 4, $\lim_{t\rightarrow\infty}N_{h}(t) = N_{h}^{*}$ for all $1\leq h\leq m$. Denote $\Phi(t, s, x_{0})$ the solution of system (3.1) with $x(s) = x_{0}$, and denote $\Theta(t, x_{0})$ the solution of system (3.2) with $y(0) = x_{0}$. Then by Proposition 1.1 in ^[18], $\Phi$ is asymptotically autonomous semiflow with limit semiflow $\Theta$. Let $\mathcal{O}_{\Phi}(s, x) = \{\Phi(t, s, x): t\geq s\}$, $x\in V$, then $\mathcal{O}_{\Phi}(s, x)$ has compact closure in $V$ since it's bounded. Let $\omega = \omega(s, x)$ which is the $\omega$-limit set of $\mathcal{O}_{\Phi}(s, x)$. By Theorem 1.8 in ^[18], we conclude that $\omega$ is non-empty, compact and connected, and it attracts $\Phi(t, s, x)$. Moreover, $\omega$ is invariant for the semiflow $\Theta$ and is chain recurrent for $\Theta$.

Now suppose that $\omega\neq\{I_{i}^{*}\}$. There exists $x = (x_{1}, x_{2}, ..., x_{m})\in \omega$ such that $x\neq\{I_{i}^{*}\}$. Let $\epsilon = d(x, \{I_{i}^{*}\})$. Since $\omega$ is compact, there exists $T > 0$ such that $d(\Theta(t, x_{0}), \{I_{i}^{*}\}) < \frac{\epsilon}{2}$ for all $x_{0}\in\omega$ and $t\geq T$. By the definition of chain recurrence, there is an $(\frac{\epsilon}{2}, T)$ chain from $x$ to $x$, i.e., there is a sequence $\{x = x_{1}, x_{2}, ..., x_{n+1} = x; t_{1}, t_{2}, ..., t_{n}\}$ for $x_{i}\in\omega$ and $t_{i}\geq T$ such that $d(\Theta(t_{i}, x_{i}), x_{i+1}) < \frac{\epsilon}{2}$. Then $d(\Theta(t_{n}, x_{n}), x) < \frac{\epsilon}{2}$, which indicates that $d(x, \{I_{i}^{*}\})\leq d(\Theta(t_{n}, x_{n}), x)+d(\Theta(t_{n}, x_{n}), \{I_{i}^{*}\}) < \frac{\epsilon}{2}+\frac{\epsilon}{2} = \epsilon$, which contradicts with $\epsilon = d(x, \{I_{i}^{*}\})$. Therefore, $\omega = \{I_{i}^{*}\}$, i.e., $\{I_{i}^{*}\}$ attracts all the forward orbits going through $V$ in system (3.1).

4.   Impact of change of productive age and birth rates on infant disease rate

In this section, we assume that the birth rate of age group $i$ is given by $b_{i}(N_{i}) = p_{i}N_{i}e^{-qN}$, we'll analyze how do demographic distribution $\frac{N_{i}^{*}}{N}$ and infant disease rate at endemic equilibrium $\frac{I_{1}^{*}}{N_{1}^{*}}$ change with birth parameters $p_{i}$, $q$ and productive age $k$.

By plugging $b_{i}(N_{i}) = p_{i}N_{i}e^{-qN}$ into Eq (2.4), we have

for $1\leq h\leq m$.

for $1 < i, h < m$. Combine the above two formulas, we have

Thus

Solving for $N$ gives

From Eq (2.4), we define

for $1\leq h\leq m$. It follows that

By plugging Eq (4.2) into $N = \sum_{h = 1}^{m}N_{h}^{*}$, we have

Equations (4.2) and (4.3) give

It's obvious from Eq (4.4) that change of $p_{i}$, $q$ and $k$ don't make a change on the demographic distribution $\frac{N_{h}^{*}}{N}$ for $1\leq h\leq m$.

Proposition 10. Assume that the birth rate of age group $i$ is given by $b_{i}(N_{i}) = p_{i}N_{i}e^{-qN}$. The stablized demographic distribution $N_{h}^{*}/N$ doesn't change with birth parameters $p_{i}$, $q$ and reproductive age $k$.

Note that change of $p_{i}$, $q$ and $k$ does have an impact on the total number $N$ and $N_{h}^{*}$, though they don't influence the ratio $N_{h}^{*}/N$.

It can be seen from Eq (4.1) that $N$ increases as $p_{i}$ increases, decreases as $q$ or $k$ increases, so does $N_{h}^{*}$ for all $1\leq h\leq m$.

Now we study the impact of changes of $p_{i}$, $q$ and $k$ on the infant disease rate $\frac{I_{1}^{*}}{N_{1}^{*}}$.

In the last section, we get that when $R_{0} > 1$, there is a endemic equilibrium of system 3.2 where $I_{1}^{*}$ satisfies $\sum_{i = 1}^{m}\frac{\lambda_{i}\frac{N_{i}^{*}}{N}N}{l_{i}+\frac{l_{1}}{M_{1}}} = 1$ where $l_{i} = \frac{\sigma+\mu_{i}}{\alpha_{i}}$ and $M_{i} = \frac{N_{i}^{*}}{I_{i}^{*}}-1$. By taking derivative with respect to $q$, we get

So

which means that $\frac{N_{1}^{*}}{I_{1}^{*}}$ increases as $q$ increases, so $\frac{I_{1}^{*}}{N_{1}^{*}}$ decreases as $q$ increases.

Similarly, we have

which means that $\frac{N_{1}^{*}}{I_{1}^{*}}$ decreases as $p_{i}$ increases, so $\frac{I_{1}^{*}}{N_{1}^{*}}$ increases as $p_{i}$ increases.

If $k$ gets larger, $N$ gets smaller as discussed above, $\frac{N_{i}^{*}}{N}$ doesn't change, since we have $\sum_{i = 1}^{m}\frac{\lambda_{i}\frac{N_{i}^{*}}{N}N}{l_{i}+\frac{l_{1}}{M_{1}}} = 1$, $\frac{l_{1}}{M_{1}}$ decreases thus $M_{1} = \frac{N_{1}^{*}}{I_{1}^{*}}$ increases. So $\frac{I_{1}^{*}}{N_{1}^{*}}$ gets smaller. By the same argument, if $k$ gets smaller, $\frac{I_{1}^{*}}{N_{1}^{*}}$ gets larger.

In conclusion, we have

Proposition 11. Assume that the birth rate of age group $i$ is given by $b_{i}(N_{i}) = p_{i}N_{i}e^{-qN}$. With all the other parameters fixed, the infant disease rate at endemic equilibrium $\frac{I_{1}^{*}}{N_{1}^{*}}$ increases as birth rate $p_{i}$ increases, and decreases as the productive age $k$ or $q$ increases.

5.   Impact of change of productive age and birth rates on disease distribution

In this section, we assume that $R_{0} > 1$ and the birth rate of age group $i$ is given by $b_{i}(N_{i}) = p_{i}N_{i}e^{-qN}$, we'll analyze how does disease distribution $\frac{I_{i}^{*}}{I^{*}}$ at endemic equilibrium change with birth parameters $p_{i}$ and $q$.

From Eq (3.6) we define $H: = \frac{\alpha_{h}(N_{h}^{*}/I_{h}^{*}-1)}{\sigma+\mu_{h}}$ for any $1\leq h\leq m$, then

and

Therefore,

Let $Q_{h}: = \sum_{j = 1\;}^{m}\frac{\alpha_{j}(\sigma+\mu_{j})}{((\sigma+\mu_{j})H+\alpha_{j})^{2}}\frac{N_{j}^{*}}{N}-\frac{\sigma+\mu_{h}}{(\sigma+\mu_{h})H+\alpha_{h}}\sum_{j = 1\;}^{m}\frac{\alpha_{j}}{(\sigma+\mu_{j})H+\alpha_{j}}\frac{N_{j}^{*}}{N}$

Proposition 12. Assume that $R_{0} > 1$ and the birth rate of age group $i$ is given by $b_{i}(N_{i}) = p_{i}N_{i}e^{-qN}$. With all the other parameters fixed, how $\frac{I_{h}^{*}}{I^{*}}$ changes with $q$, $p_{i}$ and $k$ depends on the sign of $Q_{h}$: $\frac{I_{h}^{*}}{I^{*}}$ increases as $q$ increases if $Q_{h} > 0$ and decreases as $q$ increases if $Q_{h} < 0$; $\frac{I_{h}^{*}}{I^{*}}$ decreases as $p_{i}$ increases if $Q_{h} > 0$ and increases as $p_{i}$ increases if $Q_{h} < 0$. In particular, if $1\leq h\leq m$ is such that $\frac{\alpha_{h}}{\sigma+\alpha_{h}} < \frac{\alpha_{j}}{\sigma+\alpha_{j}}$ for all $j\neq h$, $\frac{I_{h}^{*}}{I^{*}}$ decreases as $q$ increases and increases as $p_{i}$ increases; if $1\leq h\leq m$ is such that $\frac{\alpha_{h}}{\sigma+\alpha_{h}} > \frac{\alpha_{j}}{\sigma+\alpha_{j}}$ for all $j\neq h$, $\frac{I_{h}^{*}}{I^{*}}$ increases as $q$ increases and decreases as $p_{i}$ increases.

Proof. From Eq (5.3) and the conclusion we get that $\frac{N_{j}^{*}}{N}$ doesn't change with $q$, $p_{i}$ or $k$

From last section $d(\frac{N_{h}^{*}}{I_{h}^{*}})/dq > 0$, which implies that $\frac{dH}{dq} > 0$. So we have $d(\frac{I_{h}^{*}}{I^{*}})/dq > 0$ if $Q_{h} > 0$ and $d(\frac{I_{h}^{*}}{I^{*}})/dq < 0$ if $Q_{h} < 0$, which means that $\frac{I_{h}^{*}}{I^{*}}$ increases as $q$ increases if $Q_{h} > 0$ and decreases as $q$ increases if $Q_{h} < 0$. In particular, if $1\leq h\leq m$ is such that $\frac{\alpha_{h}}{\sigma+\alpha_{h}} < \frac{\alpha_{j}}{\sigma+\alpha_{j}}$ for all $j\neq h$, then $\frac{\sigma+\mu_{h}}{(\sigma+\mu_{h})H+\alpha_{h}} > \frac{\sigma+\mu_{j}}{(\sigma+\mu_{j})H+\alpha_{j}}$, which implies that $Q_{h} < 0$ thus $d(\frac{I_{h}^{*}}{I^{*}})/dq < 0$, $\frac{I_{h}^{*}}{I^{*}}$ decreases as $q$ increases; if $1\leq h\leq m$ is such that $\frac{\alpha_{h}}{\sigma+\alpha_{h}} > \frac{\alpha_{j}}{\sigma+\alpha_{j}}$ for all $j\neq h$, then $\frac{\sigma+\mu_{h}}{(\sigma+\mu_{h})H+\alpha_{h}} < \frac{\sigma+\mu_{j}}{(\sigma+\mu_{j})H+\alpha_{j}}$, which implies that $Q_{h} > 0$ thus $d(\frac{I_{h}^{*}}{I^{*}})/dq > 0$, $\frac{I_{h}^{*}}{I^{*}}$ increases as $q$ increases.

Similarly, we have

From last section we have $d(\frac{N_{h}^{*}}{I_{h}^{*}})/dp_{i} < 0$, which implies that $\frac{dH}{dp_{i}} < 0$. So we have $d(\frac{I_{h}^{*}}{I^{*}})/dp_{i} < 0$ if $Q_{h} > 0$ and $d(\frac{I_{h}^{*}}{I^{*}})/dp_{i} > 0$ if $Q_{h} < 0$, which means that $\frac{I_{h}^{*}}{I^{*}}$ decreases as $p_{i}$ increases if $Q_{h} > 0$ and increases as $p_{i}$ increases if $Q_{h} < 0$. In particular, if $1\leq h\leq m$ is such that $\frac{\alpha_{h}}{\sigma+\alpha_{h}} < \frac{\alpha_{j}}{\sigma+\alpha_{j}}$ for all $j\neq h$, then $\frac{\sigma+\mu_{h}}{(\sigma+\mu_{h})H+\alpha_{h}} > \frac{\sigma+\mu_{j}}{(\sigma+\mu_{j})H+\alpha_{j}}$, which implies that $Q_{h} < 0$ thus $d(\frac{I_{h}^{*}}{I^{*}})/dp_{i} > 0$, $\frac{I_{h}^{*}}{I^{*}}$ increases as $p_{i}$ increases; if $1\leq h\leq m$ is such that $\frac{\alpha_{h}}{\sigma+\alpha_{h}} > \frac{\alpha_{j}}{\sigma+\alpha_{j}}$ for all $j\neq h$, then $\frac{\sigma+\mu_{h}}{(\sigma+\mu_{h})H+\alpha_{h}} < \frac{\sigma+\mu_{j}}{(\sigma+\mu_{j})H+\alpha_{j}}$, which implies that $Q_{h} > 0$ thus $d(\frac{I_{h}^{*}}{I^{*}})/dp_{i} < 0$, $\frac{I_{h}^{*}}{I^{*}}$ decreases as $p_{i}$ increases.

6.   Family planning strategies

In this section, we study how family planning strategies influence demographic distribution at equilibrium, basic reproduction number and infant disease rate.

For simplicity, we assume that there are only two productive age groups, the $k$th and $(k+1)$th group. We also assume that the birth function takes a more general form $b_{i}(N_{i}) = p_{i}q(N_{k}+N_{k+1})N_{i}$ for $i = k, k+1$, where $q$ is a decreasing function. Since the maximal children each female has per unit time in age group $i$, given by $p_{i}$, are dependent on each other, more precisely, $p_{k}+p_{k+1} = b$ for some constant $b$, which is the maximal children each female has per unit time, we assume that $p_{k} = b\alpha$ and $p_{k+1} = b(1-\alpha)$, then $\alpha$ indicates the tendency to have children at an earlier age. We study how $\alpha$ influence demographic distribution at equilibrium, basic reproduction number and infant disease rate.

6.1.   Impact of family planning strategy on demographic distribution at equilibrium

From Eq (2.4), we have

Further calculation gives

where

This implies that $\alpha$ doesn't influence the demographic distribution at equilibrium.

6.2.   Impact of family planning strategy on the basic reproduction number and infant disease rate

Define

for $1 < i, k < m$. By plugging $b_{i}(N_{i}) = p_{i}q(N_{k}+N_{k+1})N_{i}$ into the equation for $N_{k}^{*}$ in Eq (6.1), we obtain

Define $F(\alpha): = b\alpha+b(1-\alpha)a_{k+1, k}$, then $F'(\alpha) = b(1-a_{k+1, k})$. From Eq (6.2), we have the following conclusions: If $a_{k+1, k} > 1$, then $F'(\alpha) < 0$, $q$ is monotone increasing with respect to $\alpha$, thus $N_{k}^{*}$ is monotone decreasing with respect to $\alpha$; if $a_{k+1, k} < 1$, then $F'(\alpha) > 0$, $q$ is monotone decreasing with respect to $\alpha$, thus $N_{k}^{*}$ is monotone increasing with respect to $\alpha$; if $a_{k+1, k} = 1$, then $F'(\alpha) = 0$, $q$ doesn't change with $\alpha$, thus $N_{k}^{*}$ doesn't change with $\alpha$.

Note that the basic reproduction number is given by

So we have the following proposition.

Proposition 13. If $a_{k+1, k} > 1$, $R_{0}$ decreases as $\alpha$ increases; if $a_{k+1, k} < 1$, $R_{0}$ increases as $\alpha$ increases; if $a_{k+1, k} = 1$, $R_{0}$ doesn't change with $\alpha$.

Proposition 13 implies that if $a_{k+1, k} > 1$, the basic reproduction number decreases when more people are having children at an early age, if $a_{k+1, k} < 1$, the basic reproduction number increases when more people are having children at an early age, if $a_{k+1, k} = 1$, the basic reproduction number doesn't depend on tendency on birth age.

From previous calculation, we have

where $l_{i} = \frac{\sigma+\mu_{i}}{\alpha_{i}}$, $M_{1} = \frac{N_{1}^{*}}{I_{1}^{*}}-1$. It implies that $\frac{I_{1}^{*}}{N_{1}^{*}}$ is monotone increasing with respect to $N_{k}^{*}$. So we have the following

Proposition 14. If $a_{k+1, k} > 1$, $\frac{I_{1}^{*}}{N_{1}^{*}}$ decreases as $\alpha$ increases; if $a_{k+1, k} < 1$, $\frac{I_{1}^{*}}{N_{1}^{*}}$ increases as $\alpha$ increases; if $a_{k+1, k} = 1$, $\frac{I_{1}^{*}}{N_{1}^{*}}$ doesn't change with $\alpha$.

Proposition 14 implies that if $a_{k+1, k} > 1$, the infant disease rate decreases when more people are having children at an early age, if $a_{k+1, k} < 1$, the infant disease rate increases when more people are having children at an early age, if $a_{k+1, k} = 1$, the infant disease rate doesn't depend on tendency of birth age.

7.   Conclusions

We proposed a stage-structured model of childhood infectious disease transmission dynamics. The population demographics dynamics is governed by a certain family and population planning strategy which gives rise to nonlinear feedback delayed effects on the reproduction ageing and rate.

The long-term aging-profile of the population is described by the pattern and stability of equilibrium of the demographic model. For this demographic model, conditions on the birth functions and death rate were given to guarantee the existence and stability of the positive equilibrium. This implies conditions on birth function and age dependent death rate to reach a stable population. We also investigate the disease transmission dynamics, using the epidemic model when the population reaches the positive equilibrium (limiting equation). We establish conditions for the existence, uniqueness and global stability of the disease endemic equilibrium and prove the global stability of the endemic equilibrium for the original epidemic model with varying population demographics. Birth function, age-dependent death rate, recover rate and transmission coefficients are all involved in these conditions.

The global stability of the endemic equilibrium allows us to examine the effects of reproduction ageing and rate, under different family planning strategies, on the childhood infectious disease transmission dynamics. We find that increasing birth rate increases the infant disease rate and reproduction aging decreases the infant disease rate. We also find that reproduction ageing and rate doesn't change the demographic distribution at equilibrium.

We investigate impacts of family planning strategies on demographic distribution at equilibrium, basic reproduction number for childhood disease and infant disease rate. We find the conditions under which planning to have a child at an early age helps to decrease/increase the basic reproduction number and infant diseases rate. We also examine demographic distribution, diseases reproductive number, infant disease rate and age distribution of disease.

For original contributions, the model we propose is new as it is stage structured and the growth through age stages is described by time delay leading to nonlinear feedback, the idea of studying the impact of population policy and family planning strategy on disease transmission dynamics is also novel. This model can be modified to fit specific childhood diseases for specific purposes. For example, it can be modified to study the impact of China's second-child policy on pertussis transmission dynamics by incorporating more compartments to distinguish children from one-child and two-children families. The work can also be potentially used to inform targeted age group for optimal vaccine booster programs.

Conflict of interest

The authors declare no conflict of interest.