A mathematical model for frogeye leaf spot epidemics in soybean

1.   Introduction

Frogeye leaf spot (FLS), caused by the fungus Cercospora sojina, is a common soybean disease ^[1,2]. First reported in Japan in 1915, the disease has spread to most soybean-growing countries throughout the world. In the United States, FLS has historically been most common in the southern region, especially the states of Alabama, Arkansas, Kentucky, Louisiana, Mississippi, and Tennessee, which account for a large portion of soybean production in the nation. Recent years have seen FLS reported with increasing frequency and severity in the midwestern and northern regions of the United States, including states such as Illinois, Iowa, Minnesota, Nebraska, North Dakota, and Wisconsin ^[3,4,5].

FLS represents a significant threat to soybean production by reducing photosynthetic leaf areas, causing premature defoliation and reduced seed weight, and eventually leading to yield losses. It is estimated that soybean yield losses due to the FLS disease can range from 30% to 60% ^[6,7]. In Argentina alone, losses due to FLS during the 2009–2010 crop season were estimated at about 2 billion USD ^[8]. In the United States and Canada, the average soybean yield losses caused by FLS were up to 1,453,225 metric tons per year during the 10-year period 2010–2019; furthermore, the estimated yield losses from 2015 to 2019 were more than twice that reported in the previous 5-year period from 2010 to 2014 ^[9,10]. Hence, effective prevention and intervention of FLS are tremendously important for agriculture, food, and the economy.

The causative agent of FLS, Cercospora sojina, can survive and overwinter in crop residue left in the soil that provides the initial inoculum ^[2,11]. Infectious conidia (i.e., asexual spores) are produced on infested residue in spring and dispersed to nearby plants via wind or rain splashes, which represents the primary route of the infection. The pathogen causes small, irregularly circular leaf lesions, typically seen in the upper canopy of the infected plants. Individual lesions can merge to form larger patterns of blight on the leaves. Such symptoms can be observed 8–12 days after infection, and new spores are produced within these lesions and subsequently spread to other plants, which represents the secondary route of the FLS disease ^[1]. Hence, the transmission and spread of FLS are polycyclic and involve multiple pathways, including both the environment (i.e., the soil) and the infected plants ^[12]. Meanwhile, some conidia generated and dispersed from infected plants will deposit on the soil surface, which may increase the primary inoculum and may contribute to the fungal growth in the soil ^[13]. This represents the feedback from the infected plants to the environment. Additionally, soybean seeds infected with Cercospora sojina may also spread the fungus when they are planted.

FLS is but one example of infectious plant diseases, which can be caused by a wide variety of pathogens including bacteria, fungi, and viruses ^[14]. Mathematical models for plant infections, although not as extensively used as for human and animal infections, started several decades ago ^[15]. Since then significant progress has been made in this area, as reviewed in ^{[16,17,18,19,20,21,22]}. Common mathematical tools for plant pathology and epidemiology include, among others: (1) Disease progress curves, which employ curves of prescribed shapes (such as monomolecular, exponential, logistic, and Gompertz functions) to describe the epidemic progression over time ^[23,24]. (2) Area under the disease progress curve, which can be used as a measure of epidemic development and may be further applied to hypothesis testing and regression analysis ^[22,25]. (3) Disease cycle models, which utilize different variables for different stages of the disease cycle, and which employ a prescribed function (with parameters fitted from data) at each stage to model the infection process ^[26,27,28]. (4) Linked differential equations, which are formulated to describe the rate of change for each state variable, typically referred to as a compartment, and the flow of infection through different compartments ^[21,29].

Among these modeling techniques, linked differential equations provide a general and powerful approach to examine the transmission, spread, and progression of plant diseases. The equations can not only characterize the generation of new infections and the transfer of individuals and exchange of information between compartments but also naturally incorporate the host population growth and the pathogen dynamics. By describing the rate of change, instead of prescribing the trajectory of evolution, for each state variable, models based on differential equations offer a more flexible way than others to investigate the complex infection process that could potentially yield a better mechanistic understanding of the disease. Compared to their applications in plant epidemiology, linked differential equations have a longer history in modeling infectious diseases for humans and other animals ^[30,31]. Classical epidemic models such as the SIR (susceptible-infectious-recovered) and SEIR (susceptible-exposed-infectious-recovered) are based on this approach. On the other hand, since analytical solutions are generally impossible to find for nonlinear differential equations, dynamical system theories are typically applied to analyze the equilibrium points and stability properties, and numerical simulation is normally conducted to approximate the solution orbits so as to gain useful insights into the infection dynamics ^{[32,33,34,35,36]}.

Simple models based on the area under the disease progress curve have been utilized for frogeye leaf spot ^[37]. Other studies (e.g., ^[38]) have applied statistical regression to predict the progression of the FLS disease. To our knowledge, however, no mathematical models based on differential equations have been published for frogeye leaf spot thus far. To fill this gap, we propose an innovative modeling study for FLS in this paper using linked differential equations. Our goal is to achieve a deep mechanistic understanding of FLS dynamics, taking into account both the disease transmission among soybean plants and the intrinsic growth and decay of the fungus Cercospora sojina in the soil. To that end, we will modify the standard SEIR model and combine it with an additional equation that represents the environmental dynamics of the pathogen, forming a coupled differential-equation system that incorporates the multiple transmission pathways of FLS. Although such a modeling technique has been utilized for some human and animal infections that involve both direct and indirect transmission routes (e.g., cholera and brucellosis) ^[39,40,41], its application to plant diseases appears to be new. We will conduct detailed equilibrium and stability analyses for the proposed model, and will use numerical simulation to verify the analytical results and to implement the model with real data for a practical application.

The remainder of this paper is organized as follows. The formulation of our mathematical model is presented in Section 2. Details of the equilibrium and stability analyses are provided in Section 3, followed by numerical simulation and data fitting in Section 4. Conclusions are drawn and some discussion is made in Section 5.

2.   Model formulation

Let $S$, $E$, $I$ and $R$ denote the densities of the susceptible (healthy), exposed, infectious, and removed (post-infectious) soybean plants, respectively, in terms of the number of plants per unit area in a field. Let also $B$ denote the density of the fungus Cercospora sojina in the soil. We propose the following model to describe the transmission and spread of frogeye leaf spot in soybean plants:

Here $N = S+E+I+R$ is the total plant density in a given field, which is assumed to be a constant. The parameter $\mu$ is the natural growth and removal rate of the plants, $\beta$ and $\alpha$ are the primary and secondary transmission rates, respectively, $\lambda$ is the reciprocal of the mean latent period, $\delta$ is the removal rate of infectious plants, $r$ is the fungal intrinsic growth rate, $k$ is the carrying capacity of the fungus, $\tau$ is the removal rate of the fungus from the soil, and $\xi$ is the average rate of contribution from an infectious plant to the fungus in the soil. All these parameters are assumed to be non-negative constants.

System (2.1) consists of nonlinear and coupled differential equations for five state variables. In this model, healthy plants (represented by the $S$ compartment) contract the infection either through the primary inoculum from the fungus in the soil, or through the secondary inoculum from the fungus produced by the infected plants, and then enter the exposed ($E$) compartment. We represent each of these two transmission modes using a bilinear incidence based on the law of mass action, where the incidence rate is directly proportional to the product of the densities of the susceptible plants and the infectious plants (or, the environmental pathogen). Exposed plants go through a latent period of the length $\lambda^{-1}$, after which they are capable of spreading the disease and enter the infectious ($I$) compartment. The last equation in system (2.1) describes the intrinsic growth of Cercospora sojina in the soil by a logistic model with the growth rate $r$ and carrying capacity $k$. Additionally, the term $\xi I$ represents the total of the reciprocal contribution per unit time from the infected plants to the environmental pathogen.

3.   Mathematical analysis

3.1.   Basic reproduction number

We conduct a detailed analysis of the disease risk, equilibrium points, and stability properties for our FLS model. First, it can be easily observed from system (2.1) that all the solutions will remain non-negative as long as the initial conditions are non-negative. Next, from the last equation in system (2.1), we have

It follows that $B$ has an upper bound $B_{max}$. Thus, we consider system (2.1) in the following positively invariant set

as a biologically meaningful domain. Clearly, there is a unique disease-free equilibrium (DFE) at

To determine the basic reproduction number ${\mathcal R}_0$ for the model, we consider $E$, $I$, and $B$ as the infectious compartments. Using the next-generation matrix technique ^[42], we obtain the nonnegative matrix $F$ representing the generation of new infections and the non-singular matrix $V$ representing the transfer of individuals between compartments:

Then we compute $V^{-1}$ and, consequently, the next-generation matrix $FV^{-1}$, as follows,

Thus, the basic reproduction number can be determined by the spectral radius of $FV^{-1}$; that is, $\mathcal{R}_0 = \rho(FV^{-1})$. Since

where $\mathbb{I}$ denotes the corresponding identity matrix, we have

Equation (3.7) shows that the disease risk of FLS, quantified by the basic reproduction number $\mathcal{R}_0$, is shaped collectively by the primary infection (represented by the terms associated with $\beta$), the secondary infection (represented by the terms associated with $\alpha$), and the environmental dynamics of the fungus (represented by the terms associated with $r$ and $\xi$).

3.2.   Nontrivial equilibrium

A nontrivial equilibrium $X = (S, E, I, R, B)$ for system (2.1) satisfies

At a nontrivial equilibrium, at least one of the variables $E$, $I$, $R$ and $B$ should be greater than $0$. Since all the variables are non-negative (to make biological sense), it is easy to observe from Eqs (3.8)–(3.12) that as long as one of these four equilibrium components is positive, all the other equilibrium components also become positive. Let $\theta = \frac{\mu N\lambda}{(\mu+\lambda)(\mu+\delta)}$. Canceling $S$ and $E$ from (3.8)–(3.10), we obtain

Solving the above quadratic equation for nonnegative $I$, we have $I = g_1(B)$, where

In addition, we have $B = p(I)$ from (3.13), where

Hence $B = p(g_1(B))$, $B\ge0$. Then we have $g_1'(B) = 1/p'(g_1(B)) = 1/p'(I)$ for $B > 0$. For an endemic equilibrium, $p(I) = B$ is positive, which implies that

It follows from (3.12) that $I = g_2(B)$, where

To obtain an endemic equilibrium, we need to solve the equation

Let $m = \min(\theta, \frac{\mu}{\alpha})$. Then $0 < \theta-I < m$ from (3.16), and $\frac{\mu}{\alpha}\theta\ge m^2$. Hence we have

Thus,

and

On the other hand,

Notice that $g_1(0) = \max\{0, \theta-\frac{\mu}{\alpha}\}\ge0 = g_2(0)$. By considering the intersection of the two curves $I = g_1(B)$ and $I = g_2(B)$ for $B > 0$, we have

(1) if $\theta > \frac{\mu}{\alpha}$, i.e., $\frac{\alpha N\lambda}{(\mu+\lambda)(\mu+\delta)} > 1$, then $g_1(0) > g_2(0)$ and Eq (3.18) has a unique solution;

(2) if $\theta = \frac{\mu}{\alpha}$, i.e., $\frac{\alpha N\lambda}{(\mu+\lambda)(\mu+\delta)} = 1$, then $g_1(0) = g_2(0)$, $g_1'(0+) = +\infty > g_2'(0)$, and Eq (3.18) has a unique solution;

(3) if $\theta < \frac{\mu}{\alpha}$, i.e., $\frac{\alpha N\lambda}{(\mu+\lambda)(\mu+\delta)} < 1$, then $g_1(0) = g_2(0)$, and

Note that

Thus,

(a) if $\mathcal{R}_0 > 1$, then $g_1'(0) > g_2'(0)$ and Eq (3.18) has a unique solution;

(b) if $\mathcal{R}_0\le1$, then $g_1'(0)\le g_2'(0)$ and there is no solution for Eq (3.18).

From Eq (3.7), $\mathcal{R}_0 > \max\left(\frac{\alpha N\lambda}{(\mu+\lambda)(\mu+\delta)}, \, \frac{r}{\tau}\right)$. It follows that $\frac{\alpha N\lambda}{(\mu+\lambda)(\mu+\delta)}\ge1$ implies $\mathcal{R}_0 > 1$. Hence, there is a unique endemic equilibrium (EE)

for system (2.1) when $\mathcal{R}_0 > 1$. The above results can be summarized as the following theorem.

Theorem 3.1. If $\mathcal{R}_0 \le 1$, system (2.1) has a unique equilibrium, the DFE. If $\mathcal{R}_0 > 1$, system (2.1) has two equilibria, the DFE and the unique endemic equilibrium (EE).

3.3.   Stability property

We establish the following stability theorems that characterize the main dynamical behavior of our FLS model.

Theorem 3.2. If $\mathcal{R}_0\le 1$, the DFE $X_0$ of system (2.1) is globally asymptotically stable in the domain $\Gamma$. If $\mathcal{R}_0 > 1$, the DFE is unstable.

Proof. Define the vector $Z = [E, \, I, \, B]^T$ for the infected compartments. One can verify that

where the matrices $F$ and $V$ are given in Eq (3.4). Introduce the vector

It follows from the fact $\mathcal{R}_0 = \rho(FV^{-1}) = \rho(V^{-1}F)$ and direct calculation that $\mathbf{u}$ is a left eigenvector associated with the eigenvalue $\mathcal{R}_0$ of the matrix $V^{-1}F$; i.e., $\mathbf{u} V^{-1}F = \mathcal{R}_0 \mathbf{u}$. Let us consider the Lyapunov function

Differentiating $L$ along the solution of system (2.1), we have

Case 1: If $\mathcal{R}_0 < 1$, then the equality $\frac{dL}{dt} = 0$ implies that $\mathbf{u} Z = 0$, which leads to $I = B = 0$ since $\mathcal{R}_0 > \frac{\alpha N\lambda}{(\mu+\lambda)(\mu+\delta)}$. The trajectory that starts in $\{(S, E, I, R, B)\in\Gamma : I = B = 0\}$ and remains in it for all $t > 0$ can only be the DFE $X_0$. That is, the largest positive invariant set on $\{(S, E, I, R, B)\in\Gamma : \frac{dL}{dt} = 0\}$ is the singleton $\{X_0\}$.

Case 2: If $\mathcal{R}_0 = 1$, i.e.,

then $1 > \max\left(\frac{\alpha N\lambda}{(\mu+\lambda)(\mu+\delta)}, \, \frac{r}{\tau}\right)$. Thus we have $(\mu+\lambda)(\mu+\delta)-\alpha N\lambda > 0$ and $\tau-r > 0$. The equality $\frac{dL}{dt} = 0$ implies that

Hence, $S = N$ and $B = 0$. Similarly, we can obtain that the largest positive invariant set on $\{(S, E, I, R, B)\in\Gamma : \frac{dL}{dt} = 0\}$ is the singleton $\{X_0\}$.

Therefore, in either case, the largest invariant set on which $\frac{dL}{dt} = 0$ consists of only the singleton $X_0 = (N, \, 0, \, 0, \, 0, \, 0)$. By LaSalle's Invariance Principle ^[43], the DFE $X_0$ is globally asymptotically stable in $\Gamma$ when $\mathcal{R}_0\le 1$.

In contrast, if $\mathcal{R}_0 > 1$, then it follows from the continuity of vector fields that $\frac{dL}{dt} > 0$ in a neighborhood of the DFE in $\mathring{\Gamma}$. Thus, the DFE is unstable based on the Lyapunov stability theory.

Before we proceed, we make two remarks here. First, based on Theorem 4.3 in ^[44], it follows that our system (2.1) is uniformly persistent when $\mathcal{R}_0 > 1$. This can be established by using the standard arguments that the DFE $X_0$ is the only equilibrium on the boundary of the domain $\Gamma$ and that $X_0$ is unstable when $\mathcal{R}_0 > 1$ (see Theorem 2.2 of ^[45] or Proposition 3.3 of ^[46]). Second, the local asymptotic stability of the endemic equilibrium $X_1$ when $\mathcal{R}_0 > 1$ can be proved by using the Routh-Hurwitz criterion, and the details are provided in Appendix A.

In what follows, we will construct a Lyapunov function to establish a stronger result and show the global asymptotic stability of the unique endemic equilibrium $X_1$ in $\mathring{\Gamma}$, the interior of $\Gamma$.

Theorem 3.3. If $\mathcal{R}_0 > 1$, the EE $X_1$ is globally asymptotically stable in $\mathring{\Gamma}$.

Proof. Let

where the symbol $Y$ can be replaced by any of the state variables $S, E, I,$ and $B$. Then $f(Y)\ge0$ for $Y, \, Y_1 > 0$ and

Now we consider a Lyapunov function in the following form:

Then $L_1\ge0$ in $\mathring{\Gamma}$ and

By using Eqs (3.8)–(3.12) and the inequality $1-x\le-\ln x$ for $x > 0$, we have

Hence,

and $\frac{dL_1}{dt} = 0$ if and only if $(S, E, I, B) = (S_1, E_1, I_1, B_1)$. Thus, the largest invariant set in $\mathring{\Gamma}$ such that $\frac{dL_1}{dt} = 0$ is the singleton $\{X_1 = (S_1, E_1, I_1, R_1, B_1)\}$. Additionally, the uniform persistence of system (2.1) ensures that for all solution orbits starting in $\mathring{\Gamma}$, we can find a compact set inside the domain such that LaSalle's Invariance Principle ^[43] can be applied. Therefore, $X_1$ is globally asymptotically stable in $\mathring{\Gamma}$.

Theorems 3.1–3.3 describe the essential dynamics of the FLS model (2.1). These results show that the condition $\mathcal{R}_0 = 1$ is a sharp threshold for the stability, where a forward transcritical bifurcation takes place. As long as $\mathcal{R}_0 > 1$, the disease will persist in the hosts (i.e., the soybean plants). In the next section, we will use numerical simulation to verify these analytical predictions and to demonstrate a real-world application of our model through data fitting.

4.   Numerical simulation

In practical field studies, the infection data for FLS are typically recorded by the rating of disease severity, ranging between $0$ to $100 \%$, which measures the percentage of currently infected soybean plants. For convenience of simulation and data fitting, we introduce the following new variables that scale the original variables $S$, $E$ and $I$ with respect to $N$, the total plant density in the given field:

We can then re-write system (2.1) as

where we have dropped the equation for $R(t)$ (the removed plants) since it is not needed in the numerical simulation. Under system (4.2), each of the variables $s(t)$, $e(t)$ and $i(t)$ ranges between 0 and 1, and the disease-free state corresponds to $s = 1$ and $e = i = B = 0$. All the simulation results presented in this section are generated from Matlab.

4.1.   Threshold dynamics

We first numerically illustrate the stability properties of the DFE and EE. Figure 1 shows a phase portrait of $i(t)$ vs. $s(t)$ when $\mathcal{R}_0 = 0.92$. Each (red) curve represents a solution orbit that is determined by a unique initial condition, which differs for different orbits. We observe that all these trajectories eventually converge to the disease-free equilibrium corresponding to $s_0 = 1$ and $i_0 = 0$, indicating the global asymptotic stability of the DFE. The pattern is similar for other scenarios with $\mathcal{R}_0 < 1$. These results confirm the analytical prediction in Theorem 3.2.

Meanwhile, Figure 2 displays a phase portrait of $i(t)$ vs. $s(t)$ for $\mathcal{R}_0 = 1.48$, where we observe that all the solution orbits converge to the endemic equilibrium at $(s_1, \, i_1)$ with $0 < s_1 < 1$ and $0 < i_1 < 1$. Other simulation results with varied values of $\mathcal{R}_0 > 1$ show a similar pattern and are not presented here. These results demonstrate the global asymptotic stability of the EE which is consistent with the prediction in Theorem 3.3.

4.2.   Data fitting and application

We now fit our model (4.2) to the disease severity data reported from ^[37] as a real-world application of our study. The field research conducted in ^[37] was designed such that each plot consisted of four rows spaced 76.2 cm apart with each row being 6 m long. Soybeans were planted at a rate of 12 per meter along each row, giving a plant density of $N \doteq 21.0$ per square meter for each plot. Soybeans were planted in May and harvested in September/October, with varied dates in different years. We take the average length of a soybean growing season as 150 days, which yields a natural growth and removal rate of $\mu = 1/150$ per day. FLS symptoms typically appear 8-12 days after a plant is infected, and we take the average latent period as $\lambda^{-1} = 10$ days. The infection can span the soybean reproductive growth stages from flowering (R1) to beginning maturity (R7), and we take the average infectious period as $\delta^{-1} = 75$ days. The fungus Cercospora sojina can survive in crop residue left in the soil for 2 years ^[2,11], which leads to a removal rate of $\tau = 1/(2*365)$ per day from the soil. Since Cercospora sojina spreads the infection through the conidia, we may interpret $B$ as the number density of the conidia in the soil environment. Based on the experimental study performed in ^[11], the mean number of Cercospora sojina conidia in the month of May was about 21,250 per gram of leaf tissue collected from the overwintering crop residue. Given that the leaf mass density normally ranges between 0.1–0.5 grams per ml and that the crop residue is relatively dry, we take 0.2 g/ml as the average leaf mass density. This yields a number density of approximately 4000 conidia per ml in May, which will be used as the initial condition for $B(t)$ in our model. Additionally, the primary inoculum concentration can be as high as 60,000 conidia per ml ^[11], which will be used as the carrying capacity $k$ in our model. The values of all these parameters are listed in Table 1.

There are, however, no published sources for the FLS primary transmission rate ($\beta$), the FLS secondary transmission rate ($\alpha$), and the rate of contribution from infectious plants to the fungus in the soil ($\xi$). These three parameters play critical roles in shaping the epidemics of FLS. We will estimate these parameters ($\alpha$, $\beta$ and $\xi$) by fitting the active infections $i(t)$ in our model system (4.2) to the FLS disease data reported in ^[37], where the disease severity was rated and recorded for each growth stage from R1 (flowering) to R7 (beginning maturity).

Before we present the fitting results, we briefly discuss the identifiability of these parameters. We apply the scale invariance local structural identifiability method, described in ^[47] and also summarized in Appendix B, to our system (4.2), where $\alpha$, $\beta$ and $\xi$ are the parameters to be analyzed, and where only the variable $i(t)$ is observed (from the FLS disease severity data). Following this method, we re-scale the non-observed variables $s$, $e$, $B$ and the parameters $\alpha$, $\beta$, $\xi$ as follows

where $u_s$, $u_e$, $u_B$, $u_\alpha$, $u_\beta$ and $u_\xi$ are the unknown scaling factors. Let $f_x$ denote the right-hand side of the equation for the variable $x$ in system (4.2), where $x = s, \, e, \, i, \, B$. We then write each $f_x$ as a summation of linearly independent functions given below:

where $f_{xj}$ denotes the $j$th summand for the function $f_x$, with $j = 1, \, 2$ and $x = s, \, e, \, i, \, B$. The invariance equations (i.e., Equation (B.5)) for our system then become

It is easy to see that system (4.5) admits only one solution

It follows that the parameters $\alpha$, $\beta$ and $\xi$ are all identifiable when the variable $i(t)$ is observed. We note, however, that this result is concerned with structural identifiability and does not necessarily indicate practical identifiability of these parameters. For more detailed discussion of parameter identifiability, we refer readers to the review article ^[48].

We list the dataset used for our model fitting in Table 2. Based on our discussion before and the data in Table 2, the initial conditions for the model (4.2) are set as $s(0) = 97\%, \; e(0) = 0, \; i(0) = 3 \%$, and $B(0) = 4000/{\rm ml}$.

Figure 3 shows the fitting result for $i(t)$ against the disease severity data presented in Table 2. The fitted values for the three parameters $\alpha$, $\beta$ and $\xi$ and their confidence intervals are listed in Table 3. Based on these parameter values, we are able to evaluate the basic reproduction number for FLS using Eq (3.7), and we obtain $\mathcal{R}_0 \approx 18.3$. The value of $\mathcal{R}_0$ for FLS appears to be high, compared to most human infections whose basic reproduction numbers typically range between 1 and 5 ^[32]. This large $\mathcal{R}_0$, fitted from field data, stems from the rapid spread of the infection and the high level of disease prevalence associated with FLS – 36% of the total host population (i.e, the soybean plants) become infected in a period of less than 5 months. The result indicates the high infection risk for FLS and the importance of effective control for this disease.

Various disease control measures for FLS have been implemented in order to improve the soybean yields ^[1]. For example, cultural practices such as tillage and crop rotation may reduce the primary infection rate $\beta$, and chemical control such as the use of fungicides may reduce the secondary infection rate $\alpha$ and the environmental fungal contribution rate $\xi$ from infected plants. To examine the effects of these intervention strategies, and to quantify the impact of these three critical parameters on model outcomes, we have conducted a series of numerical simulations with varied parameter values.

The fitted values for $\alpha$, $\beta$ and $\xi$ listed in Table 3 are regarded as the baseline values for these parameters. In the first numerical experiment, we reduce the primary infection rate $\beta$ to $50\%$, $20\%$, and $0\%$, respectively, in reference to its base value, while keeping the other parameters fixed. We conduct the simulation for each of these scenarios and plot the curves for $i(t)$ in Figure 4. We clearly observe the reduction of the disease severity with decreased values of $\beta$. In particular, when $\beta$ is reduced to $20\%$ of its base value, the disease severity is kept under $20\%$ throughout the soybean growing season. We also perform the same experiments for the secondary infection rate $\alpha$ and the fungal contribution rate $\xi$, with the results presented in Figures 5 and 6, respectively. We see a similar pattern that reducing each of these parameters to $20\%$ of its base value would push the disease severity curve under $20\%$ for the entire period. Among these, the reduction of $\alpha$ appears to have the most significant impact on lowering the disease severity (compare Figure 5 with Figures 4 and 6), indicating the importance of reducing the secondary infection from infected plants through control measures such as fungicides.

On the other hand, Figures 4–6 clearly show that changing a single parameter (among $\alpha$, $\beta$ and $\xi$) would not be sufficient to control the FLS prevalence – even with a parameter reduced to 0 so that the transmission mode associated with that parameter is totally removed, the disease severity still increases and eventually reaches near or above $10\%$. These results underscore the complexity of the multiple transmission pathways involved in FLS, indicating that using a single intervention method may not be sufficiently effective to control the FLS epidemic. Hence, we have conducted several additional sets of simulations where we assume that multiple FLS control strategies are implemented so that these parameters may be simultaneously reduced.

Figure 7 displays the results when various combinations of the parameters $\beta$, $\alpha$ and $\xi$ are simultaneously reduced to $50\%$, $20\%$, and $10\%$ of their respective base values. For each case, the improvement in disease control can be easily observed, in comparison with the single intervention results shown in Figures 4–6. The best performance is achieved when all the three parameters are simultaneously reduced (see Figure 7(d)). For cases (b), (c) and (d) in Figure 7, when the parameters are reduced to $20\%$ of their base values, the disease severity decreases from the very beginning and throughout the entire season, indicating that no FLS epidemic occurs. Changing these parameters to $10\%$ of their base values would further push the disease severity curve downward, though not significantly different from the $20\%$ scenario.

5.   Discussion

Although epidemic models based on differential equations and dynamical systems have been extensively used for infectious diseases of humans and other animals, their applications to plant diseases are less popular thus far. The mathematical model proposed in this paper for frogeye leaf spot represents a new contribution to the theoretical study of this disease and to the field of plant pathology and epidemiology. Utilizing coupled nonlinear differential equations, our model is able to describe the complex interplay between the primary infection, the secondary infection, and the environmental dynamics of the fungus, an advantage not shared by other modeling approaches such as the area under the disease progress curve and the regression analysis that have been utilized for this disease ^[37,38]. The proposed model allows us to conduct a deep investigation into the transmission, spread, and progression of FLS among soybean plants.

Mathematical analysis of this model has established a sharp threshold at $\mathcal{R}_0 = 1$ that separates two distinct types of dynamical behavior: disease eradication when $\mathcal{R}_0 < 1$ and disease persistence when $\mathcal{R}_0 > 1$. This standard result, which is applicable to a wide variety of epidemic models, indicates that prevention and intervention methods need to be implemented such that the basic reproduction number may be reduced below unity in order to eliminate the disease. We have verified these analytical predictions for our FLS model using numerical simulation.

Our data fitting and simulation studies show that the model outcomes well represent the reported disease severity data and that our model can be utilized for realistic investigation of the transmission and spread of FLS, including the impact of various disease control measures. In particular, our numerical results demonstrate that a combination of multiple intervention methods may be more effective than a single method in controlling the FLS epidemics. On the other hand, we note that the quality of our data fitting may be impacted by the limited amount of time series disease data for FLS. The accuracy of our model output could be further improved by the availability of FLS disease data with higher resolutions in the future.

For frogeye leaf spot, there are several disease control strategies that can be put into agricultural practices. These include tillage, crop rotation, resistant cultivars, seed treatment, and use of fungicides ^[1,2,37]. Tillage can help bury infested residue and decrease the chance of fungal transmission from the environment to the plants, thus reducing the primary infection. Foliar fungicides can kill the fungus generated from infected plants and prevent the growth of the spores, thus reducing the secondary infection and the reciprocal feedback from the infected plants to the environment. In addition, planting resistant soybean varieties could prevent the spread of the disease, rotation with crops not susceptible to FLS would allow time for the inoculum in the field to degrade before the next soybean planting season, and seed treatment may reduce the risk of introducing infected seeds into a field. From the modeling perspective, these intervention methods can effectively reduce the basic reproduction number associated with FLS. We have only conducted a coarse-grained simulation study to represent some of these control methods. Our FLS model can be extended to incorporate more detailed effects of these disease control measures, based on which an optimal control study can be performed to explore the most effective and practically feasible intervention strategy for FLS. We plan to pursue this direction in our future efforts.

Acknowledgments

The authors thank the four anonymous reviewers for their comments which have improved the quality of the original manuscript. The work of J.W. was partially supported by the National Science Foundation under Grant Numbers 1951345 and 2324691.

Use of AI tools declaration

The authors declare that they have not used any Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix A.   Local stability of the endemic equilibrium

Theorem A.1. If $\mathcal{R}_0 > 1$, the EE of system (2.1) is locally asymptotically stable.

Proof. Linearizing system (2.1) at the endemic equilibrium $X_1 = (S_1, \, E_1, \, I_1, \, R_1, \, B_1)$ and using Eqs (3.8)–(3.12) yield the Jacobian matrix at $X_1$:

It is easy to see that the characteristic polynomial of $J_1$ is

where $\mathbb{I}$ represents the corresponding identity matrix and

are all positive. To prove the local asymptotical stability of $X_1$, it suffices to show $(a_3a_2-a_1)a_1 > a_3^2a_0$ according to the Routh-Hurwitz criterion. For simplicity, we denote $m = 2\mu+\lambda+\delta$, $M = m+\frac{\mu N}{S_1}+\frac{rB_1}k$, and $\theta = \frac{\beta\lambda\mu B_1N}{I_1}+\alpha\lambda\mu(N-S_1)$. We may express $a_3^2a_0$ as a function of $\xi$ in the following way,

where

Similarly, we can also express $(a_3a_2-a_1)a_1$ as a function of $\xi$,

We note that $\theta > \lambda\mu(N-S_1)(\frac{\beta B_1}{I_1}+\alpha)$ and

After some tedious algebraic calculations, we find

Thus, we have $(a_3a_2-a_1)a_1 > a_3^2a_0$ for any $\xi\ge0$, as expected. This completes the proof for the local asymptotic stability of the endemic equilibrium $X_1$.

Appendix B.   The scale invariance local structural identifiability method

In what follows, we summarize the scale invariance local structural identifiability method introduced in ^[47]. Consider a system of differential equations involving $n$ variables $x_j$ ($j = 1, \cdots, n$) and $m$ parameters $\lambda_k$ ($k = 1, \cdots, m$):

It is assumed that the variables $x_1, \, \cdots, \, x_r$ can be observed (measured) from the experiment, while the variables $x_{r+1}, \, \cdots, \, x_n$ cannot be observed.

Write each function $f_j$ as a summation of linearly independent terms:

where $\widetilde{x}_k$ and $\widetilde{\lambda}_k$ denote the subset of variables and parameters, respectively, included in the function $f_{jk}$ for $1 \leq k \leq M$ and $j = 1, \cdots, n$. The method consists of the following steps:

1) Scale the unobserved variables and parameters by

where $u_{x_j}$ and $u_{\lambda_k}$ are unknown scaling factors associated with the variable $x_j$ and the parameter $\lambda_k$, respectively. Substitute the scaled variables and parameters into the equations

for $j = 1, \cdots, n$.

1) Equate each linearly independent function $f_{jk}$ to its scaled counterpart; i.e.,

where $u_{x_j} = 1$ for $1 \leq j \leq r$. These are referred to as the invariance equations.

3) Solve the invariance equations for the scaling factors $u_{x_j}$ and $u_{\lambda_k}$ or find combinations of these scaling factors that leave the system invariant.

4) Any parameter $\lambda_k$ associated with a solution $u_{\lambda_k} = 1$ is identifiable. Any variable $x_j$ associated with a solution $u_{x_j} = 1$ is observable. In contrast, those parameters whose scaling factors are coupled will form identifiable groups but cannot be identified independently.