Model formulation
Using pigs in the large-scale farm as nodes in the network and considering the type of pig and pig house, age, and corresponding transfer, an individual-based state probability model is established to study the effectiveness of immune programs. The time unit of the model is taken as day.
The set of pig types in the farm is denoted by Z={c,r,s,b}, where c represents commercial pigs, r represents reserve pigs, s represents breeding sows, and b represents breeding boars. We sort and number all pigs, pig houses and pens. The serial number set of all pigs in the farm is I={1,2,3,...,n}, where n is the total number of pigs. The serial number set of all pig houses in the farm is H={1,2,3,...,m}, where m is the total number of pig houses; the serial number set of the delivery room is H1={1,2,3,...m1}; the serial number set of the nursery house is H2={m1+1,m1+2,m1+3,...m2}; the serial number set of the fattening house is H3={m2+1,m2+2,m2+3,...m3}; the serial number set of the reserve house is H4={m3+1,m3+2,m3+3,...m4}; the serial number set of the pregnancy room is H5={m4+1,m4+2,m4+3,...m5}; and the serial number set of the boar breeding station is H6={m5+1,m5+2,m5+3,...m}. The serial number set of all pens in pig houses is F={1,2,3,...,g}, where g is the total number of pens; the serial number set of pig pens in the delivery room is F1={1,2,3,...,g1}; the serial number set of pig pens in the nursery house is F2={g1+1,g1+2,g1+3,...,g2}; the serial number set of pig pens in the fattening house is F3={g2+1,g2+2,g2+3,...,g3}; the serial number set of pig pens in the reserve house is F4={g3+1,g3+2,g3+3,...,g4}; the serial number set of pig pens in the pregnancy room is F5={g4+1,g4+2,g4+3,...,g5}; and the serial number set of pig pens in the boar breeding station is F6={g5+1,g5+2,g5+3,...,g}. For the i-th (i∈I) pig, namely, the i-th node, zi(t) denotes the type of the i-th pig at time t, and zi(t)∈Z; ai(t) denotes the days old of the i-th pig at time t, and ai(t)≥0; hi(t) denotes the serial number of the pig house where the i-th pig is located at time t, and hi(t)∈H; and fi(t) denotes the serial number of the pen where the i-th pig is located at time t, and fi(t)∈F.
When pigs are infected with FMD by direct contact or environmental transmission, pigs need 3−4 weeks to clear the virus. During the infection period, some infected pigs are culled due to clinical symptoms. The remaining pigs can be self-healing. According to this transmission mechanism, each pig may experience four states, including susceptible, immunized, infected and recovered, which are expressed as S, V, E and R, respectively. Therefore, for the i-th node, the probabilities of being in five states at time t are \(P_{S}^{(z_{i}(t),a_{i}(t),h_{i}(t),f_{i}(t))}(i,t)\), \(P_{V}^{(z_{i}(t),a_{i}(t),h_{i}(t),f_{i}(t))}(i,t)\), \(P_{E}^{(z_{i}(t),a_{i}(t),h_{i}(t),f_{i}(t))}(i,t)\), and \(P_{R}^{(z_{i}(t),a_{i}(t),h_{i}(t),f_{i}(t))}(i,t)\), respectively. In addition, the birth, transfer, and sale of pigs are reflected in the superscript (zi(t),ai(t),hi(t),fi(t)). The above processes of birth, sale, transfer of pigs and state transfer due to the transmission of FMD can be simulated by Python software. At time t, the amount of FMDV in the environment of the l-th pig house is represented by B(l,t), where l∈H.
The variation of the above 4 probabilities with the course of infection and immunization and the state probability transfer diagram of the i-th pig are shown in Fig. 7, and system (1) is its corresponding system.
$$ \left\{\begin{array}{l}{a}_i(t)={a}_i\left(t-1\right)+1,\\ {}\begin{array}{cc}{P}_S^{\left({z}_i(t),{a}_i(t),{h}_i(t),{f}_i(t)\right)}\left(i,t\right)\kern0.3em =& k(t){P}_V^{\left({z}_i\left(t\kern0.3em -\kern0.3em 1\right),{a}_i\left(t\kern0.3em -\kern0.3em 1\right),{h}_i\left(t\kern0.3em -\kern0.3em 1\right),{f}_i\left(t\kern0.3em -\kern0.3em 1\right)\right)}\left(i,t\kern0.3em -\kern0.3em 1\right)\kern0.3em +\kern0.3em \Big(1-{\lambda}_1^{f_i\left(t\kern0.3em -\kern0.3em 1\right)}\left(i,t\right)\\ {}\kern0.3em -\kern0.3em {\lambda}_2^{h_i\left(t\kern0.3em -\kern0.3em 1\right)}\left(i,t\right)\kern0.3em -\kern0.3em \delta (t)\Big){P}_S^{\left({z}_i\left(t\kern0.3em -\kern0.3em 1\right),{a}_i\left(t\kern0.3em -\kern0.3em 1\right),{h}_i\left(t\kern0.3em -\kern0.3em 1\right),{f}_i\left(t\kern0.3em -\kern0.3em 1\right)\right)}\left(i,t\kern0.3em -\kern0.3em 1\right),\\ {}{P}_V^{\left({z}_i(t),{a}_i(t),{h}_i(t),{f}_i(t)\right)}\left(i,t\right)\kern0.3em =& \left(1-k(t)\right){P}_V^{\left({z}_i\left(t\kern0.3em -\kern0.3em 1\right),{a}_i\left(t\kern0.3em -\kern0.3em 1\right),{h}_i\left(t\kern0.3em -\kern0.3em 1\right),{f}_i\left(t\kern0.3em -\kern0.3em 1\right)\right)}\left(i,t\kern0.3em -\kern0.3em 1\right)\\ {}+\kern0.3em \delta (t){P}_S^{\left({z}_i\left(t\kern0.3em -\kern0.3em 1\right),{a}_i\left(t\kern0.3em -\kern0.3em 1\right),{h}_i\left(t\kern0.3em -\kern0.3em 1\right),{f}_i\left(t\kern0.3em -\kern0.3em 1\right)\right)}\left(i,t\kern0.3em -\kern0.3em 1\right),\\ {}{P}_E^{\left({z}_i(t),{a}_i(t),{h}_i(t),{f}_i(t)\right)}\left(i,t\right)\kern0.3em =& \left(1\kern0.3em -\kern0.3em \eta \right){P}_E^{\left({z}_i\left(t\kern0.3em -\kern0.3em 1\right),{a}_i\left(t\kern0.3em -\kern0.3em 1\right),{h}_i\left(t\kern0.3em -\kern0.3em 1\right),{f}_i\left(t\kern0.3em -\kern0.3em 1\right)\right)}\left(i,t\kern0.3em -\kern0.3em 1\right)\kern0.60em +\kern0.3em \Big({\lambda}_1^{f_i\left(t\kern0.3em -\kern0.3em 1\right)}\left(i,t\right)\\ {}\kern0.3em +\kern0.3em {\lambda}_2^{h_i\left(t\kern0.3em -\kern0.3em 1\right)}\left(i,t\right)\Big){P}_S^{\left({z}_i\left(t\kern0.3em -\kern0.3em 1\right),{a}_i\left(t\kern0.3em -\kern0.3em 1\right),{h}_i\left(t\kern0.3em -\kern0.3em 1\right),{f}_i\left(t\kern0.3em -\kern0.3em 1\right)\right)}\left(i,t\kern0.3em -\kern0.3em 1\right),\\ {}{P}_R^{\left({z}_i(t),{a}_i(t),{h}_i(t),{f}_i(t)\right)}\left(i,t\right)\kern0.3em =& \upmu \upeta {P}_E^{\left({z}_i\left(t\kern0.3em -\kern0.3em 1\right),{a}_i\left(t\kern0.3em -\kern0.3em 1\right),{h}_i\left(t\kern0.3em -\kern0.3em 1\right),{f}_i\left(t\kern0.3em -\kern0.3em 1\right)\right)}\left(i,t\kern0.3em -\kern0.3em 1\right)\\ {}+{P}_R^{\left({z}_i\left(t\kern0.3em -\kern0.3em 1\right),{a}_i\left(t\kern0.3em -\kern0.3em 1\right),{h}_i\left(t\kern0.3em -\kern0.3em 1\right),{f}_i\left(t\kern0.3em -\kern0.3em 1\right)\right)}\left(i,t\kern0.3em -\kern0.3em 1\right),\\ {}\end{array}\\ {}B\left(l,t\right)=\left(1-\omega -\nu \right)B\left(l,t-1\right)\kern0.3em +\kern0.3em \sum \limits_{h_j\left(t\kern0.3em -\kern0.3em 1\right)\kern0.3em =\kern0.3em l}^{j\in I}\alpha {P}_E^{\left({z}_j\left(t\kern0.3em -\kern0.3em 1\right),{a}_j\left(t\kern0.3em -\kern0.3em 1\right),{h}_j\left(t\kern0.3em -\kern0.3em 1\right),{f}_j\left(t\kern0.3em -\kern0.3em 1\right)\right)}\left(j,t\kern0.3em -\kern0.3em 1\right).\end{array}\right. $$
(1)
In this model, two modes of transmission, direct contact transmission and environmental transmission, are considered. The probabilities of the i-th pig being infected by the two methods at time t are expressed as \(\lambda _{1}^{f_{i}(t-1)}(i,t)\) and \(\lambda _{2}^{h_{i}(t-1)}(i,t)\), respectively. The specific expressions for \(\lambda _{1}^{f_{i}(t-1)}(i,t)\) and \(\lambda _{2}^{h_{i}(t-1)}(i,t)\) and meanings of other parameters in system (1) are interpreted as follows.
(1) \(\lambda _{1}^{f_{i}(t-1)}(i,t)\). Direct contact can only occur in the same pen. The probability of a susceptible pig being infected by an infected pig in the same pen through direct contact per unit time is denoted by β. Therefore, the probability that the i-th pig cannot be infected by the j-th pig in the same pen per unit time is
$$1\,-\,\beta P_{E}^{(z_{j}(t\,-\,1),a_{j}(t\,-\,1),h_{j}(t\,-\,1),f_{j}(t\,-\,1))}(j,t\,-\,1); $$
The probability that the i-th pig is not infected by all pigs in the same pen is
$${\prod_{\substack{j\in I, j\neq i \\ f_{j}(t-1)=f_{i}(t-1)}}\left[1-\beta P_{E}^{(z_{j}(t-1),a_{j}(t-1),h_{j}(t-1),f_{j}(t-1))}(j,t-1)\right].} $$
Therefore, the probability of the i-th pig being infected through direct contact is
$$\begin{array}{*{20}l} {}\lambda_{1}^{f_{i}(t-1)}(i,t)&=\!1\!-\prod_{\substack{j\in I, j\neq i \\ f_{j}(t-1)=f_{i}(t-1)}}\\ &{}\quad\big[\!1\,-\,\beta P_{E}^{(z_{j}(t-1),a_{j}(t-1),h_{j}(t-1),f_{j}(t-1))}(j,t\,-\,1\!)\big]. \end{array} $$
(2) \(\lambda _{2}^{h_{i}(t-1)}(i,t)\). For environmental transmission, it is necessary to consider not only the transmission caused by media or aerosols within pig houses but also the transmission caused by media and human activities among pig houses. Here, γ is the infection rate coefficient of virus in the environment to one susceptible pig, and p is the probability that pigs are exposed to viruses in the environment in different pig houses. At time t, the probability that the i-th pig is infected by the virus in the environment in the same pig house is γB(hi(t − 1),t − 1); the probability that the i-th pig is infected by the virus in the environment in different pig houses is \(\sum _{\bar {h}\neq h_{i}(t\,-\,1)}^{\bar {h}\in H} p\gamma B(\bar {h},t\,-\,1)\). Therefore, at time t, the probability that the i-th pig is infected by FMDV in the environment is
$${\lambda_{2}^{h_{i}(t-1)}(i,t)=\gamma B(h_{i}(t-1),t-1)+\!\sum_{\substack{\bar{h}\in H \\ \bar{h}\neq h_{i}(t-1)}}p\gamma B(\bar{h},t-1).} $$
(3) k(t): Antibody decay rate.
(4) δ(t): Antibody growth rate.
(5) η: Virus clearance rate.
(6) μ: The probability of self-recovery.
(7) α: The amount of virus excreted by the infected pig per unit time.
(8) ω: The natural decay rate of the virus in the environment per unit time.
(9) ν: The disinfection rate of the virus in the environment per unit time.
If the infection process of the farm is not considered and only immunization, production and sale are considered, the system (1) becomes the following: [!t]
$$\begin{array}{@{}rcl@{}} \left\{\begin{array}{ll} a_{i}(t)=a_{i}(t-1)+1,\\ P_{S}^{(z_{i}(t),a_{i}(t),h_{i}(t),f_{i}(t))}(i,t)\!=&k(t)P_{V}^{(z_{i}(t\,-\,1),a_{i}(t\,-\,1),h_{i}(t\,-\,1),f_{i}(t\,-\,1))}(i,t\,-\,1)\!\\ &+\!(1-\!\delta(t))P_{S}^{(z_{i}(t\,-\,1),a_{i}(t\,-\,1),h_{i}(t\,-\,1),f_{i}(t\,-\,1))}(i,t\,-\,1),\\ P_{V}^{(z_{i}(t),a_{i}(t),h_{i}(t),f_{i}(t))}(i,t)\!=&(1-k(t))P_{V}^{(z_{i}(t\,-\,1),a_{i}(t\,-\,1),h_{i}(t\,-\,1),f_{i}(t\,-\,1))}(i,t\,-\,1)\\ &+\!\delta(t)P_{S}^{(z_{i}(t\,-\,1),a_{i}(t\,-\,1),h_{i}(t\,-\,1),f_{i}(t\,-\,1))}(i,t\,-\,1).\\ \end{array}\right. \end{array} $$
(2)
At time t, the number of pigs whose immune antibody is qualified in the whole farm is \(\sum _{i\in I}P_{V}^{(z_{i}(t),a_{i}(t),h_{i}(t),f_{i}(t))}(i,t)\), the immune antibody qualified rate in population is \(\frac {\sum _{i\in I}P_{V}^{(z_{i}(t),a_{i}(t),h_{i}(t),f_{i}(t))}(i,t)}{N(t)}\), the positive number of pigs being infected is \(\sum _{i\in I}P_{E}^{(z_{i}(t),a_{i}(t),h_{i}(t),f_{i}(t))}(i,t)\), and the positive rate in population is \(\frac {\sum _{i\in I}P_{E}^{(z_{i}(t),a_{i}(t),h_{i}(t),f_{i}(t))}(i,t)}{N(t)}\).
Parameter values and initial values
The values and sources of parameters in system (1) are interpreted as follows:
(1) Pigs infected with FMD require 3−4 weeks to remove the virus (Alexandersen et al. 2002). So, \(\eta =\frac {1}{28}\).
(2) According to the actual situation, no pigs with symptoms were culled at the farm from 2016 to 2019. So, μ=1.
(3) The average peak amount of virus discharged by pigs per day can reach 106.1TCID50 (Zhang et al. 2019). So, α=106.1.
(4) The survival time of the virus in the environment is approximately 30 days (Zhang et al. 2019). So, \(\omega =\frac {1}{30}\).
(5) The farm is generally sterilized once a month. Each time disinfection rate is assumed to reach 100%. Thus, \(\nu =\frac {1}{30}\).
(6) In Additional file 3, we establish an ordinary differential dynamical system (system (0.1)) for the farm to describe the spread of FMD within the pig population. By applying the data fitting of system (0.1) with real data, β=0.9948 is obtained, and its biological meaning is the product of the probability that a susceptible pig is infected by an infected pig per contact and the number of contacts in a year. The biological meaning of β in the system (1) is the probability that a susceptible pig is infected by direct contact with an infected pig in the same pen within one day. Given that the average number of pigs per pen is 15, the average number of direct contacts is 15. Thus, for system (1), β=0.9948÷365÷15≈0.0002.
(7) In Additional file 3, the time unit of the ordinary differential dynamic system is year, and γ=4.5029×10−10. The time unit of system (1) is day. Therefore, for system (1), γ=4.5029×10−10÷365≈1.23×10−12.
(8) Assume that p=0.25.
(9) According to Additional file 4, maternal antibodies, primary immunization and secondary immunization correspond to different k(t) and δ(t). The maternal antibody of piglets reached a maximum at 7 days old. The antibody level of the first and second immunizations reached a maximum at 28 days after vaccine injection. When the immune antibody protection rate increases, k(t)=0. When the immune antibody protection rate decreases, δ(t)=0. Commercial pigs are vaccinated at T1 and T2 days of age. Therefore, according to the fitting results of the antibody fluctuation rule in Additional file 4, when zi(t)=c,
$$\begin{array}{*{20}l} k(t)&=\left\{ \begin{array}{lcl} 0, && {1\leq a_{i}(t)\leq 7,}\\ 0.0154, && {7< a_{i}(t)\leq T_{1},}\\ 0, && {T_{1}< a_{i}(t)\leq T_{1}+28,}\\ 0.0212, && {T_{1}+28< a_{i}(t)\leq T_{2},}\\ 0, && {T_{2}< a_{i}(t)\leq T_{2}+28,}\\ 0.0043, && {a_{i}(t)> T_{2}+28,} \end{array} \right.\\ &\text{and}\\ \delta(t)&=\left\{ \begin{array}{lcl} 0.2299, & & {1\leq a_{i}(t)\leq 7,}\\ 0, & & {7< a_{i}(t)\leq T_{1},}\\ 0.027, & & {T_{1}< a_{i}(t)\leq T_{1}+28,}\\ 0, & & {T_{1}+28< a_{i}(t)\leq T_{2},}\\ 0.0727, & & {T_{2}< a_{i}(t)\leq T_{2}+28,}\\ 0, & & {a_{i}(t)> T_{2}+28.} \end{array} \right. \end{array} $$
It is assumed that the increase and decrease in the antibody protection rate in the reserve pig, breeding boar and breeding sow after vaccination is the same as that in piglets after secondary immunization. Reserve pigs are vaccinated at T3 and T4 days of age. Therefore, when zi(t)=r,
$$\begin{array}{*{20}l} k(t)&=\left\{ \begin{array}{lcl} 0.0043, & & {a_{i}(t)\leq T_{3},}\\ 0, & & {T_{3}< a_{i}(t)\leq T_{3}+28,}\\ 0.0043, & & {T_{3}+28< a_{i}(t)\leq T_{4},}\\ 0, & & {T_{4}< a_{i}(t)\leq T_{4}+28,}\\ 0.0043, & & {a_{i}(t)> T_{4}+28,} \end{array} \right.\\ &\text{and}\\ \delta(t)&=\left\{ \begin{array}{lcl} 0, & & {a_{i}(t)\leq T_{3},}\\ 0.0727, & & {T_{3}< a_{i}(t)\leq T_{3}+28,}\\ 0, & & {T_{3}+28< a_{i}(t)\leq T_{4},}\\ 0.0727, & & {T_{4}< a_{i}(t)\leq T_{4}+28,}\\ 0, & & {a_{i}(t)> T_{4}+28.} \end{array} \right. \end{array} $$
Breeding boars and breeding sows are vaccinated thrice a year, and the set of immunization times is Tv={January 20, May 20, September 20}. Therefore, when zi(t)∈{s,b},
$$k(t)=\left\{ \begin{array}{lcl} 0, & & \text{within\ 28\ days \ after\ vaccination,}\\ 0.0043, & & \text{other\ time,} \end{array} \right. $$
and
$$\delta(t)=\left\{ \begin{array}{lcl} 0.0727, & & \text{within\ 28\ days \ after\ vaccination,}\\ 0, & & \text{other\ time.}\\ \end{array} \right. $$
In addition, according to the actual situation of the farm considered in this paper, T1=60, T2=90, T3=169 and T4=259.
(10) By applying system (1), the data in 2019 are used as the initial value to predict the qualified rate of immune antibodies and the level of pathogen positive rate in the pig population within one year. According to the survey in 2019, the immunization antibody qualification rate in the pig population was 89.16%. The pathogen positive rate was 0.0083; thus, the proportion of infected pigs (E) was 0.83%. Assume that the number of pigs that recover (R) is 0 and the number of viruses (B) in the environment is 0. According to Table S2 in Additional file 1, the total number of pigs in the farm (N) is 15063. It can be calculated that there were 1508 susceptible pigs, 13430 immune pigs and 125 infected pigs in the farm.
(11) For system (2), it is assumed that all pigs in the farm at the initial moment are susceptible. Thus, the probability of any pig being susceptible is 1, and the probability of being immunized is 0. We simulate the system (2) to consider the influence of different primary immunization times, secondary immunization times and immunization frequencies on the immune antibody qualified rate in the farm. The immune antibody level of breeding sows directly determines the maternal antibody level of piglets.
The control reproduction number
For one pig in a certain type of pig house, its individual basic reproduction number (Keeling and Grenfell 2000) is
$${}R_{0}=\int_{0}^{\infty}[a_{1}\beta+a_{2}\gamma h(B)+(N-a_{2})p\gamma h(B)]T{f}(T)\mathrm{d}T, $$
and its individual control reproduction number is
$${}R_{c}\,=\,\int_{0}^{\infty}\![a_{1}\beta p_{1}+a_{2}\gamma h(B)p_{1}+(N-a_{2})p\gamma h(B)p_{2}]\!T{f}(T)\mathrm{d}T. $$
where (1) a1 is the number of pigs in the same pen. p1 is the proportion of susceptible pigs in the same pig house. (2) Let the right end of \(\frac {dB}{dt}=\alpha E-(\omega +\nu)B\) in system (0.1) in Additional file 3 be equal to 0; then, \(B=\frac {\alpha E}{\omega +\nu }\). Therefore, at steady state, the amount of virus released into the air by a pig is \(h(B)=\frac {\alpha }{\omega +\nu }\). a2 is the number of pigs in the same house. N−a2 is the number of pigs that are not in the same pen. p2 is the proportion of susceptible pigs in the whole farm. (3) T is the infection period and obeys an exponential distribution, namely, T∼Exp(η+d). Then, the density function is f(T)=(η+d)e−(η+d)T, and the average infection period is \(\int _{0}^{\infty }T{f}(T)\mathrm {d}T=\frac {1}{\eta +d}\). Therefore,
$$\begin{array}{*{20}l} R_{c}&=p_{1}\left(\frac{1}{\eta+d}\right)\left(\beta a_{1}+\frac{\gamma \alpha a_{2}}{\omega+\nu}\right)\\ &\quad+p_{2}\left(\frac{1}{\eta+d}\right)\left(\frac{p\gamma \alpha (N-a_{2})}{\omega+\nu}\right), \end{array} $$
which is used as an evaluation index.