A nested compartmental model to assess the efficacy of paratuberculosis control measures on U.S. dairy farms

Paratuberculosis, also known as Johne's disease (JD), is a chronic contagious disease, caused by Mycobacterium avium subsp. paratuberculosis (MAP). The disease is incurable, fatal and causes economic losses estimated to exceed 200 million dollars to the U.S. dairy industry annually. Several preventive and control measures have been recommended; however, only a few of these measures have been validated empirically. Using a nested compartmental (NC) modeling approach, the main objective of this research was to identify the best combination of control and preventive measures that minimizes the prevalence and incidence of JD and the risk of MAP occurrence in a dairy herd. The NC model employs both MAP transmission estimates and data on pen movement of cattle on a dairy to quantify the effectiveness of control and preventive measures. To obtain reasonable ranges of parameter values for between-pen movements, the NC model was fitted to the movement data of four typical California dairy farms. Using the estimated ranges of the movement parameters and those of JD from previous research, the basic reproduction number was calculated to measure the risk of MAP occurrence in each pen environment as well as the entire dairy. Although the interventions evaluated by the NC model were shown to reduce the infection, no single measure alone was capable of eradicating the infection. The numerical simulations suggest that a combination of test and cull with more frequent manure removal is the most effective method in reducing incidence, prevalence and the risk of MAP occurrence. Other control measures such as limiting calf-adult cow contacts, raising calves in a disease-free herd or colostrum management were less effective.

Introduction Johne's disease (JD) is a chronic, infectious gastrointestinal disease of domestic and wild ruminants (i.e. cattle, sheep, goats, deer, and bison), caused by Mycobacterium avium subsp. paratuberculosis (MAP). Johne's disease is a global disease, which was first observed in dairy cows in 1895 [1]. Environmental viability studies found that MAP can survive for 8 months in feces at ambient conditions [2] and for 19 months in water at 38 degrees of centigrade. MAP remains viable in a desiccated state for up to 47 months [3]. MAP is difficult to eradicate from a herd, because the pathogen persists in the environment for a long time.
Upon infection with MAP, cattle undergo an asymptomatic period which can last for years. As the disease advances, an infected cow eventually becomes overtly symptomatic with decreased milk production, persistent diarrhea, and despite having no changes in appetite, the affected animal will exhibit progressive wasting and death if not culled. The very slow progression of JD and the difficulty in identifying infected animals due to imperfect diagnostic tests contribute to the difficulty in conducting MAP control measure studies. Accurate and detailed data from testing animals throughout their lifespan and at slaughter are often not available [4]. As a result, researchers have relied on mathematical and statistical models to study transmission dynamics of MAP. Infectious disease modeling and simulation can also be used to determine the effects of control policies and to identify the risk factors contributing to disease spread [5,6,7,8].
Although between-pen cattle movement is an important factor in the spread of infectious diseases, it has been largely ignored in the modeling and analysis of various infectious diseases of farmed ruminants. Many studies [7,8] include a diffusion term to capture the cattle movement, however, these models assume random movement of cattle throughout the farm rather than the purposeful movement of cattle between pens for management reasons. In addition, recent studies have used network data to study the dynamic of cattle trade movements [9]. However, most of network models only investigate the mobility patterns of individual animals among farms and ignore between-pen cattle movements within each farm. Without considering between-pen cattle movements it would be difficult to provide meaningful inferences pertaining to within and between-pen disease transmission dynamics.
The objective of this study was to create a nested compartmental (NC) model for MAP transmission which accounts for progression of the disease and also the movement of cattle between pens on a dairy. The model fitted to the cattle movement data and employed to assess control and preventive measures. Using the numerical simulations of the model we identified the most optimal intervention combinations yielding the greatest reduction of JD incidence and prevalence on a dairy.

NC model justification and design
Cattle movements between pens can modify the contact rates between MAP in the environment, susceptible and infected cattle. The changes in cattle population and contact patterns are more complicated than concept of random mixing assumed in most infectious disease models. In practice, the continuous management-based changes in pen populations on a dairy can strongly influence MAP transmission dynamics and the effects of MAP control measures. Thus, an NC model for MAP transmission on a dairy farm was constructed by considering the states of MAP infection in dairy cattle as well as the frequency and patterns of cattle movements between pens on a dairy. The Cattle Movement (CM) model is a compartmental model formulated with a set of Ordinary Differential Equations (ODEs) and acts as the shell or outer layer of the NC model. The CM model represents different pens on a dairy and within each pen type, a compartmental MAP transmission model is embedded according to the age of cattle in the respective pen.
The cattle movement model. The CM model encompasses the different production stages common to large dairy herds (!500 milking cows) which represented 47% of US dairies in 2007 [10]. The CM model was constructed by dividing the dairy into pen types that closely represent cattle housing on a dairy. The pens include the calf nursery, growing and breeding pens, post-calving (fresh), high and low milk production, and non-milking (dry) cow pens for a total of 14 pen types. Table 1 lists the pen types, descriptions and the expected residence times per given animal. Fig 1 is a flowchart of the CM model depicting the dynamics of moving dairy cattle between pens on a dairy. Although several dairy farms may have different pen structures, the illustrated flowchart is widely accepted among researchers and practitioners. Moreover, the flowchart represents the pen types, and there can be multiple pens of each type.
Calves born to nulliparous females also known as springers (female cattle which have never given birth before) in pen 5, are transported to the pre-weaned calf hutches, collectively identified here, pen 1. Also, calves born to parous females in pen 14, which are either uniparous (first lactation) or multiparous (second and greater lactation), are transported to pen 1. Once calves are weaned they are moved to pen 2, post-weaned group pens. On a dairy that raises female heifers for breeding, heifer calves are moved to pen 3 at breeding age and to pen 4 when pregnant. As the pregnant heifer approach calving (birth) they are moved to the springers pen, where they calve; or, depending on the dairy's management, may be moved to an individual or group maternity pen. Because it is common for springers to calve in the closeup pens, the CM model was designed as such.
The recently calved females, now first lactation (uniparous) dams, are then moved to pen 6 after calving. Pen 6 (fresh/hospital) houses only recently calved first lactation cows separated from multiparous cows. This is a common management practice since older cows tend to be more dominant and may limit the younger smaller uniparous females' access to feed. Furthermore, depending on the dairy's management, pen 6 may also serve as a hospital pen to facilitate segregation of colostrum and milk produced in the first few days post-calving since both are not saleable for human consumption and hence are not milked into the bulk tank on the dairy.
Pen 8 houses high producing (high-milk production) cows. When milk production begins to decline, first lactation cows are moved to the low-production pen 9. Later in lactation and as milking cows approach calving, they will undergo dry-off, an industry term referring to the voluntary cessation of milking approximately 60 days prior to calving, a period necessary to replenish a cow's body reserves and initiation of colostrum production in preparation for the new born calf. Dry-off cows are moved to pen 12 until before calving (usually 1 to 2 weeks) when they are then moved to the close-up pen (pen 13) and fed a different ration. Cows that start calving are moved to the calving (maternity) pen 14 and once calved they are moved to pen 7, the fresh/hospital pen.
Pen 7 houses multiparous (calved more than once) fresh cows (an industry term for female cattle that recently calved) to segregate and collect their colostrum and the first 1 to 7 days of milk, depending on the type of antibiotics used at dry-off to avoid antibiotic residues in milk sold for human consumption. As with uniparous females, high-milking cows are moved to pen 10 and as milk production declines they are moved to pen 11, the low milking pen. Moving cattle between pens includes redundancy in the fresh/hospital, high-milking, and lowmilking pens, essentially separating calves (pens 1 and 2), heifers (pens 3 to 5) and adult cows (pens 6 to 14). Thereafter, a cow circulates through these last six pens: 12, 13, 14, 7, 10, and 11, throughout its productive lifespan, which is about 4.8 years [11]. Moving cattle between pens was assumed to occur at a constant rate of d ij from pen i to pen j. Cows are moved from one pen to another in a prescribed sequence that depends on their age (i.e., calf (i = 1, 2), heifer (i = 3, 4, 5), and adult (i = 6,. . ., 14)), their state of productivity (high milking, low milking, dry cow (i = 10, 11, 12), and/or state of health or fertility (breeding, pregnancy, calving, fresh/hospital (i = 3, . . ., 7)). Birth and purchase rate (b j ), as well as culling and all-cause mortality rate (m i ) affect the pen totals at any point in time. Following the CM model structure (Fig 1), the corresponding set of ODEs was developed as equation (1) Estimating the between-pen movement rates. Data from four California dairies was accessed retrospectively through their dairy herd improvement software (DairyComp 305, Valley Agricultural Software, Tulare, CA) (S4 Appendix). Records from varying intervals of time between January 2011, to June 2015 were used to estimate rates of moving cows between pens. Dairy 1 contributed pen movement data from two different time periods demarcated by the dairy herd's transition from an all Holstein to a mixed breed (Holstein and Jersey) essentially acting as two dairies and hence bringing the total to five herds.
Herd managers and veterinarians were interviewed to identify the pen types and pen population demographics including age, production and reproductive status. Although 14 pen types were identified, cattle of the same age, production, or reproductive state were housed in one or more physical pens. Hence cow movement rates were estimated for pen type and not for each pen. Pen population records from the study dairies were examined by programming a Matlab code. The code employs the optimization toolbox (with program "fmincon.m") to estimate the time intervals and calculate each cow's unique residence time in a pen. To account for the possibility that a cow may have been moved into a pen before the record extraction date, 5 days were added to the beginning and end of each cow's unique record date, respectively. Subsequently, the rates of moving cows between pens were estimated as the inverse of the pen residence times for each dairy farm.
MAP transmission models. The progression of MAP infection is embedded into each movement model state and hence essentially in each pen's environment. However, different MAP transmission models were specified depending on cattle age. The mathematical models included a Susceptible-Latent-Environment (SLE) model in pre-weaned calves (i = 1), a Susceptible-Latent-Infectious-Environment (SLIE) model for cattle from weaning to the first calving (i = 2, . . ., 6), and a Susceptible-Latent-Infectious-Super shedder -Environment (SLICE) model for older cattle (i = 7, . . ., 14). Table 2 summarizes the variables of the above-mentioned models according to the rate of fecal MAP shedding rate as described by [12].
Specifically, the Susceptible (S i ) animals in pen i can be exposed to the infection by direct host-to-host contact, or indirectly by contacting the contaminated pen environment (i.e., the pathogen load in barn surfaces denoted by P i ) or the general environment (i.e., the pathogen load in the recycled lagoon water, commonly used to flush pen surfaces, denoted by E). The host-to-host contact occurs when susceptible individuals come into contact with infectious individuals including super-shedders. The exposed animals in pen i are infected, but yet to be infectious, hence are latent and denoted by L i . As infection progresses the latently infected cattle in pen i become infectious (I i ) and eventually may become super-shedders (C i ).  depicts the MAP transmission model in preweaned calves in pen 1 of the CM model, which distinguishes between three classes: Susceptible (S i ), Latent (L i ), Environment (E) here onwards referred to as the SLE model and the corresponding system of ODEs (2) is given in S3 Appendix. Intrauterine transmission from infected dams to their calves while rare, occurs even when a dam is subclinical [13]. This process has been considered in SLE model (Fig 2), where a is the proportion of the newborn calves that are born latently infected and 1-a is the remaining susceptible proportion. Although calves are presumably more susceptible to MAP infection compared to adult cattle [14], it is very difficult to identify infected calves due to the disease's prolonged latent period. Furthermore, studies suggest that MAP infected pre-weaned calves do not commonly shed the bacterium and the disease in this age group may not necessarily include an infectious state [15]. Hence, the SLE model assumes that the number of shedding calves is either negligible or the amount of shedding does not sustainably influence the transmission dynamics of JD (Fig 2). Furthermore, animals may show no clinical symptoms of disease for years after infection [16] and diagnostic tests are not sensitive enough to identify infected animals in this latent stage. Therefore, susceptible preweaned calves may progress to the latent stage but may not contaminate the pen environment (i.e., P 1 (t) = 0 for all t = 0).
For SLIE and SLICE models, the main assumption is that a susceptible host can become infected after direct contact with contaminated environment, an infectious host or a supershedder. Infected cattle shed the MAP bacilli into their feces and hence the pen environment, which contaminates the general environment; that infection also spreads due to the cattle moving dynamics, further exposing susceptible cattle on the dairy farm. The primary transmission route for the disease is fecal-oral [17]. Free-living MAP can survive more than a year in the environment [14]. In addition, MAP has been found in milk and colostrum [18], semen [19], blood and saliva [20]. The SLIE model is considered in pens i = 1, . . ., 6 with the assumption that calves and heifer MAP transmission do not include the super-shedding stage of MAP infection due to their younger age (Fig 3). The SLICE model is considered in pens i = 7, . . ., 14, which includes supper-shedders. See S3 Appendix for the set of ODEs (3) corresponding to the SLIE model. Assembling the NC model. The progression of MAP infection at various stages of infection was incorporated within each pen of the CM model. Calf population in pen 1 is divided into susceptible and latent individuals, with disease dynamics based on the SLE model. In Fig  1, the rates m 1 and b 1 represent the mean culling/all-cause mortality rate, and purchase/birth rate, respectively. A proportion a 1 of newborn calves are born latently infected and the rest (i.e., 1a 1 ) are susceptible. The coefficient β 1G represents the transmission rate due to exposure to MAP in the general environment including due to recycled lagoon water from the entire dairy and used to flush below the calf hutches, and 1/r is the duration of pathogen survival. Use of fresh water to flush is recommended, however recycled lagoon water collected from the parlor and flush from adult pens is sometimes used exposing calves to MAP from the remaining herd [21].
Next, we embedded the SLIE model into the CM model for pens i = 2, . . ., 6 (from weaning to calving). These pens do not include super-shedders (i.e., C i (t) = 0 for all t ! 0) due to age of the cattle. See the compartmental diagram in Fig 3 and the set of ODEs (3) in S3 Appendix, corresponding to the SLIE model.
Each of the rates m i , μ i , b i , for i = 2, . . ., 6, represents the mean culling and all-cause mortality rate, farm animal removal rate, and purchase/birth rate, respectively. Each of the coefficients β 2I , β 2P , and β 2G represents the transmission rate due to infectious cattle (β 2I ), pen environment (β 2P ), and the general environment due to recycled lagoon water used to flush the entire dairy (β 2G ). Each of σ 2L , σ 2I , U 2I , v 2 , and r represents infection state rates for change from the latent stage (σ 2L ), infectious (σ 2I ), mean infectious shedding (U 2I ), transition rate from a pen to the general environment (v 2 ), and the duration of pathogen survival (r).
The SLICE model was embedded in pens i = 7, . . ., 14 (adult cattle including super-shedders C i (t) see Fig 4). Dynamics of this model is described by the system of ODEs (4) in S3 Appendix. The SLICE model represents the progression of MAP infection in the late stages of the disease which includes Tiwari et al's [22] description of a fourth phase being the advanced clinical infection when cattle may show emaciation as body proteins are metabolized resulting in death due to cachexia. Cows are usually culled before the super-shedder state, aided by being fecal test positive with high concentrations of MAP shed in feces [22]. The assembled NC model includes the CM model with the embedded MAP transmission models, SLE, SLIE and SLICE models, which is presented in Fig 5. The NC model is formulated with three systems of ODEs (5)-(7) provided in S3 Appendix.
Control and prevention of MAP transmission. The NC model was simulated to investigate the effectiveness of the MAP control and prevention strategies with respect to cattle movements between the pens. The main objective was to identify the optimal combination of strategies for MAP control and prevention. Five strategies for prevention and control of MAP on dairy farms were contrasted. An explanation of these measures (control 1-5) is as follows.
1. Colostrum management: Feeding colostrum replacer (CR) vs. maternal colostrum (MC) is expected to reduce exposure to MAP either shed in colostrum or contaminating it at harvesting, transport or feeding to newborn calves [23].
2. Offsite heifer-rearing: Newborn calves are removed at birth with the goal of reducing exposure of calves to MAP [24,25]. Calves are instead raised off-site at a calf nursery with no adult cattle or recycled lagoon water and are returned to the source dairy at an older age such as post-weaning, breeding age, or as springers prior to their first calving. 3. Intensive environmental cleaning: MAP bioburden is reduced in the environment to decrease the natural challenge of susceptible cattle. Literature on this approach is scarce but such a measure is employed by more frequent scraping of fecal slurry in pen floors and power washing enclosure surfaces. Amongst examples of cleaning the environment of calves is the practice of cleaning hutches after weaning and before newborn calves are housed in them. Hutches are commonly cleaned by being soaking using a specially fabricated sprinkler system for hours to days before being power-washed, inverted for exposure to sunlight's ultraviolet rays for days and on some premises hutches are then sprayed using lime, also known as liming or white washing. An example of cleaning the environment of adult cattle include the scraping and cleaning of maternity pen surface either using power washing, application of disinfectants including lime, and repacking with bedding material between calvings. In the NC model, the amount of pathogen concentration in the environment can be expressed as a function of the cleaning efficiency [26]. As mathematically shown in [5,27] an exponential pathogen reduction such as 10-fold or more is possible for a farm that resumes/adapts a cleaning policy or uses fresh water instead of the lagoon water.
4a. Test and cull scenario (a): This measure relies on testing dairy cows at dry off using serum enzyme-linked immunosorbent assay (ELISA) on a weekly basis and culling testpositive cows and hence affects cows in pen 12 only [28]. A sensitivity of 34.2% and specificity of 95.8% has been measured for MAP ELISA [28,29]. We incorporated these values in the NC model simulations related to test and cull scenario.
4b. Test and cull scenario (b): It is based on testing all the adult cows (lactating and dry) annually using serum ELISA and hence affects cows in pens 7-14 [28] using the same ELISA diagnostic accuracy as in strategy 4a.

Delayed exposure:
The exposure of susceptible adult cattle to MAP infected herd mates is delayed. This is in contrast with control 2, where exposure of susceptible calf to MAP infected herd is prevented at birth. Estimates for such a measure were based on adult cow infection upon introduction to an infected herd [30].
The basic reproduction number R 0 . The basic reproduction number, R 0 , is defined as the mean number of secondary infections caused by a typical infected individual introduced into a totally susceptible population [6,31]. In particular, infection will gradually disappear if R 0 < 1; whereas, an outbreak is expected when R 0 > 1. Due to the long incubation period, variable severity and duration of clinical disease, the term "occurrence" was used instead of "outbreak" of MAP infection. The next-generation matrix approach [6,31,32] was used to estimate R 0 according to the NC model. Particularly, Infectious (I), General Environment (E), Pen Environment (P), and Super-shedder (C) compartments are considered as disease compartments, and the largest nonnegative eigenvalue of K matrix is the R 0 .
The SLE disease model of pre-weaned calves had the disease-free equilibrium (DFE) given Note that the SLE model does not have any endemic equilibrium and R ½SLE 0 ¼0; which implies that the DFE is stable, and there will be no occurrence of disease in pen 1 at any time. This is due to the fact that there is no infectious individual considered in pen 1 and the general environment can only make the individuals latent. Another way to show the above result is through the linear stability analysis of DFE.
For the SLIE Model of MAP transmission in post-weaned, breeding, pregnant, and springer pens, the DFE and the basic reproduction number are given by Similarly, for the SLICE Model of MAP transmission in adult pens, the DFE and the basic reproduction number are given by See S1 Appendix for the details of the calculations. For the NC model, the R 0 calculation was done numerically. Although full R 0 expression of the NC model was not achievable, it was possible to obtain R 0 expression according to each pen and the pens before and after it. These R 0 expressions are listed as follows: where i, j = 2,. . .,5;d i−1,j and d i,j+1 are the cattle movement rates from the pen before and to the pen after, respectively. The remaining R 0 expressions (pen [6][7][8][9][10][11][12][13][14] are listed in S3 Appendix. As mentioned above, correspond to SLE, SLIE, and SLICE models, respectively. In that case, for each i = 2,. . .,14, if R ½i 0 > 1, then there will an occurrence of MAP. However, if R ½i 0 < 1 for all i = 2, . . ., 14 does not guarantee that the disease will die out. Therefore between-pen cattle movements can trigger the occurrence and persistence of MAP infection on a dairy farm. Global uncertainty analysis. The control measures are the strategies to control or reduce the occurrence of MAP as reflected by the basic reproduction number of the NC model. In the NC model, R 0 has 63 parameters. To understand how R 0 is affected by these parameters, we completed a global uncertainty analysis [33]. A total of 27 combined or individual control strategies were designed and each was simulated 50,000 times to examine their effects on prevalence and incidence of MAP infection and R 0 value on a dairy farm of 10,000 susceptible cows with a super shedder and an infectious cow introduced to the farm. Specifically, using the range of parameter values in using the parameter values in Table 3, Tables B and C (see S3 Appendix), the numerical simulations of the NC model were carried out to estimate R 0 values and MAP incidence and prevalence. The simulation randomly picked numbers within the given ranges for each parameter used to calculate R 0 . For each case of control strategy, 50,000 calculations of R 0 values, minimum, maximum, mean, 95% confidence intervals (CI), and risk of MAP occurrence were calculated. The risk was calculated as the fraction of the simulation iterations, where R 0 was greater than 1. Table 3 includes the range parameter values used in the model simulations. In addition to S2 Appendix, details of parameter estimations are provided in the rest of this section. As mentioned before, we assumed that the main route of MAP transmission in the calf population (i.e., pen 1) is the through the general environment. Although there are some studies that consider MAP transmission from transient shedders to susceptible calves [7,15,33], we assumed that the number of calf-to-calf transmission in pen 1 is negligible. Furthermore, since the calf population is separated from the heifer and adult populations, the adult-to-calf transmission rates are considered zero. As suggested in [15,38,43], we assumed that a portion The estimate was assumed to be the highest annual percentage of infected heifers since it spanned a range of 12 to 21 months of follow up and for β G in the adult population we used the range considered in [8].
The infectious cattle transmission rate β I in the heifer population is adopted from [8]. However, β I in the adult population is calculated from the values provided in columns 3 and 4 of Table 2 in [30]. The value of β C in the calf and heifer population is considered zero due to the fact that the super-shedders are only considered in the adult population. The values of β P may vary considerably, depending on the quality of cleaning practices (e.g., scraping and power wash). The ranges of stage durations 1/σ L, 1/σ I and 1/σ C are mainly adopted from the previous studies [16,38,39]. The average shedding rate of infectious cattle in the heifer and adult population is obtained from our clinical study [12] and the work by Bolton et al. [37]. Duration of pathogen survival is considered the same for all pens. Most studies do not characterize the duration based on the host population. We therefore considered a range of 0.8 to 1.5 year as suggested by Magonbedze et al. [8]. Pathogen transmission rate from pen environments to the general environment may vary based on the age group due to the amount of manure produced by each age group. We assumed larger values for adult populations compared to heifer and calf populations. The animal removal rates vary from farm to farm, but it is known that the rates are higher in the calf and heifer populations. We adopted the same range of values used in [36].
Historical computerized dairy farm records with information on cattle IDs in each of the farms' pens were obtained with the herd veterinarians' and farm owners' permission and no animals were enrolled for the purpose of the current study.

Analysis of cow movement data
Dairy herd movement records exported from each study herd included cow identification numbers, date of record and pen location in Table A (see S3 Appendix), for a snapshot of the dairy farm data from 2011 to 2015. Some dairies had intervals of missing data for one or more pens, the latter could be due to the dairy management not utilizing such a pen or pens, or simply due to missing backups at regular intervals. The missing data was removed from the study. As summarized in Table 4, the cattle movement data of four dairy farms were used in this

Effectiveness of JD control measures.
In summary, serum ELISA test and cull is the most effective single control measure in reducing MAP infection. By far the best outcome is obtained by combining three control measures of test and cull, cleaning, and isolating calves and heifers from the herd. The risk of MAP occurrence was calculated by dividing the number of iterations with R 0 greater than one by total number of iterations that R 0 has been calculated. When we compare the no control option versus all combined control strategies, the risk of MAP occurrence in dairy cattle drops from 82% to 42% and the mean R 0 value drops from 3.92 to 0.89. Although this demonstrates a very effective approach to JD management on a dairy farm, it reveals that the MAP occurrence risk is not eliminated even though that all control measures are simultaneously applied.
Despite 42% risk of MAP occurrence, simulations of a MAP infected herd showed that employing all control measures reduced mean prevalence of MAP below 0.02% in calves and heifers, and mean prevalence in adult cows of 1.05% over ten years. Hence complete eradication of MAP was not possible, despite the fact that the prevalence and incidence of MAP were extremely low in the window of ten years. Table 5 summarizes results of the NC model simulations assuming that each control measure separately applied to a dairy farm. S1 Fig depicts the distributions of R 0 values. In each panel, the curve represents the fitted generalized extreme value distribution. See Table F (S3  Appendix) for the estimated sigma and mean values. In the absence of any control measures, identified as Control 0 in Table 5, the mean R 0 value was 3.92 and with a long tail in the frequency plot such that it exceeded R 0 = 20. The numerical simulations indicated that none of the controls were individually effective and hence they each failed to reduce the mean R 0 values to below 1. In this regard, the top three measures were controls 4b, 3 and 4a with the mean R 0 values of 1.31, 1.51, and 2.11, respectively. Table 6 summarizes statistics of R 0 values for MAP transmission under all possible binary combinations of the control measures. Also, S2 and S3 Figs shows the related R 0 distributions. Although the mean R 0 values was reduced from 3.92 to 3.85 but the combination of controls 1 and 2 were not successful in reducing the risk of infection and hence MAP transmission. Control measure 4a, weekly test and cull of cows at dry-off (pen 12), and 4b, annual test and cull of adult cows (pens 7-14) made up the most effective binary combination control measure, while a combination of Control 3 with Control 4a or 4b was the second most effective combination control measure. Nevertheless, none of these binary combinations reduced the mean R 0 value below one. However, as shown in Table 6 the risk significantly dropped in the cases in which test and cull (i.e., Control 4a or 4b) was combined with a control measure other than controls 1 or 2 (see S3 Appendix for more details of distribution curves and the related probability density functions). In the best case scenario, the risk of MAP occurrence decreased from 81% (control 0) to 47% (controls 4a and 4b). Also, in all cases the mean R 0 value was greater than 1, which indicates that JD will gradually spread in the herd.
Although risk of MAP infection decreased with triple and all control measures, the risk was not eliminated and remained non-zero. In particular, as shown in Table 7, the risk of MAP infection decreased from 81% (control 0) to 42% (controls 3, 4a and 4b) and the mean R 0 value decreased from 3.92 to 0.89. From the most effective to the least, the combined measures with Risk is the proportion of the number of times that R 0 was greater than 1.  Hence, on average, a farm that employs all control measures, or a combination of the three control measures should remain or gradually become disease free. However, there is still more than 40% risk of MAP infection remaining in the herd. In a similar manner, prevalence and incidence of MAP infection were estimated based on simulations for applying single control measures. Fig 6 represents the asymptotic behavior of MAP infection prevalence and incidence in a dairy farm. Namely, the curves were obtained by taking the mean values of 50,000 NC model simulations for a long period of time (t = 1,000 years). For each simulation, a super shedder and an infected cow are initially introduced to the herd. It can be seen that control 4b (i.e., annual test and cull of adult cows) was the most  successful method in the population of adult cows. Control 3 (i.e., intensive environmental cleaning) could effectively slow down the increase of incidence and prevalence in calf and heifer populations. It should also be noted that control 5, which was designed for delaying the exposure of calves and heifers to infected cows, was the most effective method in the calf and heifer population to keep both prevalence and incidence low. This is despite the fact that control 5 had a poor efficacy of 81% risk of MAP occurrence ( Table 5). Details of the asymptotic values associated with the single measures can be found in Table F (S3 Appendix). Further simulations indicated that a combination of test and cull (i.e. control 4a or 4b) with control 1 (i.e. Colostrum management) did not reduce the incidence and prevalence in the calf and heifer populations due to the fact that test and cull is rarely applied to calf and heifer pens.
We also calculated the range of the mean prevalence and incidence estimates in a practical time period (i.e., in the interval of 10 years) using 50,000 NC model simulations. These values are presented in Table 8, where control 5 is still the best control measure in the population of calves and heifers. In the population of adult cows there was no single control measure, which was superior to all other measures.
Hence, we investigated the prevalence and the incidence for the cases that more than one measure was employed. Namely, the results of combined control measures are presented in Table 9 showing that the mean prevalence and incidence estimates were substantially less in the calf and heifer populations.  (Table 9), result in prevalence estimates below 0.01% (one out of 10,000 cows), which is an important result that a disease-free herd can remain disease free, under two assumptions. First, extremely low number of infectious cow or supper shedder (n < 5) are accidently introduced to the herd.
Second, a combination of the above-mentioned control measures is strictly implemented. A common control measure among these effective combinations is control measure 5 under which calves are born and raised in uninfected herds delaying to exposure to infected cattle. Despite being a different scenario, such Estimates may serve as a conservative (worst case) scenario resulting in a prevalence of 0.008% (less than 1 in 10,000) estimated in the heifer population, making it the most effective measure in this age group. The next most effective control  Although the offsite nursery prevents contact between the calves and adult cattle, the environment in the off-site nursery pens may be contaminated by lagoon water in case of recycling flush water; and therefore, a mean of 0.07% prevalence is expected. Nevertheless, the simulations suggested that control 2 was the second best measure to reduce the JD prevalence in the calf and heifer populations. The mean incidence and prevalence values were extremely low due to the fact that the model simulations assumed that only one super-shedder adult cow and another infected adult cow were introduced into the herd in pen 10 and pen 8, respectively and followed for 10 years. Similar results were obtained when small numbers (i.e. n 5) of infected cows and supershedders were introduced into a herd of 10,000 cows. This indicates that the illustrated results are consistent for small numbers of infectious cows and supper shedders initially introduced to the herd.
In the population of adult cows, controls 2 & 4a, 2 & 4b, 3 & 4b, 5 & 4b, 4a & 4b (Table 9), and all controls combined (Table 10) result in a MAP prevalence of 0.52%. Measures 4a (weekly test and cull of dry cows) or 4b (annual test and cull of adult cattle) are common to all the adult cattle effective control measures. Hence an effective way to reduce MAP prevalence in the adult cow population is test and cull of test-positive cattle. However, control 4a was more effective than 4b (Table 8) resulting in a MAP prevalence of 0.61% and 0.98%, respectively. Table 10 shows the number of weekly incidence and the mean MAP prevalence for the most effective triple combination control measures and all of the control measures by the end

Conclusions
The simulation results indicate that no single control measure was sufficient to prevent increase in incidence of JD; however, Control 4b (i.e., test and cull of adult cattle in a dairy herd annually) resulted in the best single control measure. The most effective combination of binary control measures was produced by controls 4a annual test and cull of adult cows (pens 7-14) and 4b (i.e., weekly test and cull of cows at dry-off cows, pen 12). The overall risk of MAP occurrence was substantially reduced when test and cull was combined with intensive enclosure cleaning to reduce MAP concentration in the environment. Particularly, the best triple control measures resulted when combining Controls 3, 4a and 4b, which combined increased scraping of fecal slurry on solid surfaces in the dairy and /or power washing by 10-fold to reduce the environmental pathogen load, while also testing and culling dry-off cows on weekly basis and adult cattle annually. A farm that employs all control measures or a combination of these three control measures has the minimum risk of JD occurrence. It also has extremely prevalence and incidence provided that the number of infectious cow and supper shedder added to the herd is very small (i.e., less than 5 in a population of 10,000 cows).
Finally, it should be noted that these results can be expected if the dairy manager adheres to a cattle movement pattern between pens which maintains a degree of isolation between calves and cows and within the cow population as illustrated in the Cattle Movement diagram. Purposefully moving cattle between pens in a prescribed sequence changes the contact patterns between susceptible and infected cows beyond the assumption of random mixing inherent in infectious disease models. Cattle movement management is integral to the effectiveness of MAP control measures and changes to this system can modify the anticipated success of the control measures.

Limitations of the study
Modeling JD with effects of vaccination has been addressed in previous works (see for example [15,33]). In the present study, we did not investigate the effects of vaccination in our modeling and numerical simulations. Vaccination has its own shortcomings and is not practiced on several dairy farms. Previous research has shown that exposure to MAP vaccines or M. avium antigens can result in false positive tuberculosis tests, which is a concern for herds in TB free states and specifically those that commonly transport cattle across state lines [43,44,45,46]. Furthermore, no vaccine has been developed to fully protect calves [33]. There is currently no available approved treatment in food animals once an animal has contracted the MAP infection. For such reasons vaccinating against MAP is not widely practiced and hence was not considered in the current model. In the present work we assumed that the amount of shedding in the calf population does not sustainably influence the transmission dynamics of JD, i.e., γ C = 0 for pen 1 (Table 3). Nevertheless, this could be oversimplifying assumption in cases that the shedding rate is greater than a critical value.
Although the simulation result indicate that test and cull can be an effective control measure, there are two major concerns regarding test and cull. First, test and cull result is an immediate economic loss, which may not be recovered for a long period. Second, diagnostic tests to identify infected cows (e.g., ELISA-based JD control) often have low sensitivities and are often costly to apply routinely. Therefore, the efficacy of test and cull substantially varies based on the frequency and sensitivity of the test. There are simulating models [14,33,47] and field [16,48] studies that aim to determine the optimal culling rate in different herds based on the longterm profitability of the control measure. However, more data and model simulations are needed to develop reliable, effective and profitable JD control programs.
It should be noted that the data related to this study is from California dairies. Hence, the outcomes of current study may not necessarily apply to non-intensive dairy systems elsewhere in the US and the world. However, for dairies that manage cows in housing units and groups similar to the study dairies our findings may apply in terms of effectiveness of control measures and what may be expected in reduction of MAP transmission. Another limitation of the current study as with other mathematical modeling studies and specifically those modeling MAP transmission is the lack of precise transmission rates and other inputs needed by the model. Such model inputs require specifically designed studies that can limit variability and target the specific rate of interest. However, MAP's chronicity increases the duration of such studies which may translate to increase in cost in addition to prolonged duration of studies and potential for loss of follow up of study animals given other competing risks. To address these limitations, the current study identified several key assumptions that can be justified to utilize ranges of transmission rates from previous works (Table 4 and S2 Appendix).

Additional remarks
Investigating the optimal use of the cattle movement model with additional controls can benefit from these findings as the data shows that test and cull strategies seem to give the best outcome for R 0 . When test and cull is applied in pens 7 through 14 we see the most desirable outcome. While the primary goal of this work was to determine the efficacy of control measures using a NC model applied to JD on dairy farms, such models could also be employed to explore impacts on other animals and potentially applied to other diseases.  Table 7 (Table A), R 0 expressions for the NC model (pen [6][7][8][9][10][11][12][13][14], Estimated parameter values of the CM model (Tables B-C), Number observations and between-pen movements (Tables D-E), and Estimated sigma and mean values with 95% CI (Table F)