Study of an Epidemiological Model for Plant Virus Diseases with Periodic Coefficients

Abstract

In the present article, we research the existence of the positive periodic solutions for a mathematical model that describes the propagation dynamics of a pathogen living within a vector population over a plant population. We propose a generalized compartment model of the susceptible–infected–susceptible (SIS) type. This model is derived primarily based on four assumptions: (i) the plant population is subdivided into healthy plants, which are susceptible to virus infection, and infected plants; (ii) the vector population is categorized into non-infectious and infectious vectors; (iii) the dynamics of pathogen propagation follow the standard susceptible–infected–susceptible pattern; and (iv) the rates of pathogen propagation are time-dependent functions. The main contribution of this paper is the introduction of a sufficient condition for the existence of positive periodic solutions in the model. The proof of our main results relies on a priori estimates of system solutions and the application of coincidence degree theory. Additionally, we present some numerical examples that demonstrate the periodic behavior of the system.

1. Introduction

The mathematical modeling of plant virus diseases has been increasing in recent decades since it constitutes one of the most relevant tools for understanding virus propagation and implementing strategies of control, which include the development of appropriate technologies [1,2,3,4,5,6]. The interaction between plants and viruses is complex; consequently, mathematical models are built to provide information about the most relevant variables in the environment. Additionally, there is a range of techniques used to construct these mathematical models. These include regression analysis, mixed model and generalized linear model analyses, multivariate and time series analyses, differential equations, and stochastic models [6,7,8,9,10,11]. Specifically, in this paper, we focus on applying ordinary differential equations to model the dynamics of plant virus disease epidemics that affect plant populations and are transmitted by vector populations. For a more detailed review of plant virus disease epidemics we refer to [12,13,14,15].

Periodicity is a common characteristic in biologically relevant phenomena. Typically, these systems exhibit a range of scenarios, from disease-free states to mutual extinction, often including stable periodic orbits. These orbits can be influenced by seasonal variations, leading to fluctuating periods of high and low virus transmission rates. Understanding these patterns is crucial for predicting virus outbreaks and devising effective control measures. For instance, studies have highlighted the significance of seasonal variations in disease dynamics, underscoring the need for the strategic implementation of control measures during critical periods, such as the rainy season [16]. This concept of periodic environments is well-documented [17], particularly in plant virus contexts where roguing infected plants is a known control strategy [18].

For the agricultural industry, which plays a pivotal role in food security, understanding these dynamics is vital. Viral outbreaks in plants can severely impact agriculture. Periodic forcing studies are, therefore, instrumental for farmers and agricultural policymakers in strategizing crop rotations, planting schedules, and timing treatments effectively. The distinction between continuous and impulsive roguing control is a pertinent example [19]. A deeper comprehension of the periodic nature of plant virus interactions paves the way for more sophisticated control strategies. This study examines a mathematical model that incorporates periodic forcing, offering more precise predictions of virus spread and outbreak intensity. Such models are indispensable for preparing and responding to agricultural crises, ensuring a proactive approach to safeguarding crop health and food security.

Periodic forcing is especially important in scenarios where the plants are cassava plantations, the vectors are Bemisia tabaci whiteflies, and the virus is the cassava mosaic virus: Whiteflies inflict damage on a variety of plants—including cassava, tomatoes, sweet potatoes, roses, and houseplants—from their nymph stage, which lasts about 3 days, through to their adult phase, spanning 5 to 28 days. While mature adults typically fly short distances, usually no more than 2 meters during daylight, air currents play a crucial role in dispersing them across the cultivation area. Immediately upon reaching adulthood, these insects mate, with each female laying between 80 and 300 eggs [20,21]. The species Trialeurodes vaporariorum is more suited to temperate climates, whereas Bemisia tabasi thrives in warmer or tropical environments. Consequently, the emergence of these vector insects, as well as their interactions with healthy plants, occurs during times of the year when temperatures favor the presence of whiteflies [22]. This underscores the importance of investigating periodic forcing in the plant virus dynamic. Identifying the specific times of the year when whiteflies are most likely to appear, according to models, is crucial. This knowledge will facilitate the implementation of targeted mitigation and control strategies, reducing the need for indiscriminate use of insecticides and pesticides, which are detrimental to both the environment and the health of those applying them.

We begin by considering the basic dynamics of a population of plants and vectors, organisms that can transmit a pathogen [23]. We assume that the plant population is categorized into two groups: healthy plants and diseased plants, whose time-dependent density populations are given by S and I, respectively. Then, following the standard compartmental modeling [24,25,26,27], we deduce that the basic dynamics are described by the following system:

$$\begin{array}{ccc}\hfill \frac{dS}{dt}& =& \overline{r}S-{\overline{k}}_{1}SI,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \frac{dI}{dt}& =& -\overline{a}I+{\overline{k}}_{1}SI,\hfill \end{array}$$

(2)

where $\overline{r}$, ${\overline{k}}_1$, and $\overline{a}$ are the positive constants modeling the rates of replanting of healthy plants, the “contact” between healthy and diseased plants, and the death of diseased (or removal of) plants, respectively. However, the “contact” between healthy and diseased plants is produced by a vector population, and the term ${\overline{k}}_{1}SI$ should be improved by considering the vector population. We consider that the vector population is divided into non-infective and infective populations, with densities U and W, respectively. Then, if we consider that ${\overline{k}}_1$ is the contact between healthy plants and infective vectors, the term ${\overline{k}}_{1}SI$ in (1) and (2) must be replaced by ${\overline{k}}_{1}SW$. Moreover, we consider two kinds of behaviors in the dynamics of the vector population: logistic population growth and the acquisition of infectivity by the vectors through contact of non-infective individuals with diseased plants. Thus, the mathematical model is given by the following system:

$$\begin{array}{ccc}\hfill \frac{dS}{dt}& =& \overline{r}S-{\overline{k}}_{1}SW,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill \frac{dI}{dt}& =& -\overline{a}I+{\overline{k}}_{1}SW,\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill \frac{dU}{dt}& =& \overline{b}(U+W)\left[1-\frac{U+W}{\overline{m}(S+I)}\right]-{\overline{k}}_{2}IU-\overline{c}U,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill \frac{dW}{dt}& =& {\overline{k}}_{2}IU-\overline{c}W,\hfill \end{array}$$

(6)

where $\overline{b},\overline{c},{\overline{k}}_2$, and $\overline{m}$ are positive constants denoting the rates of birth, death, and contact between infected plants and non-infective vectors, and the maximum vector abundance, respectively.

Nowadays, there are more advanced mathematical models for plant virus diseases (see, for instance, [18,23,28,29,30,31,32,33,34,35]). Between them, we are interested in two models analyzed in [35]. We remark that the models given in [35] are focused on cassava plantations infested by begomoviruses carried by the whitefly Bemisia tabaci, commonly known as cassava mosaic virus. However, the model can be applied to other types of crops where a cutting is a living part (a stem, a branch, or a shoot) that has been taken from a plant in order to graft it onto another plant or into a container for further development. This multiplication method, using cuttings, is a common reproductive process for plants like fruit trees (blueberries, grapes, citrus, kiwi, cherry, etc.), various trees (ficus, olive trees, elms, willow, laurel, etc.), and ornamentals (roses, magnolias, lavender, lilacs, etc.). To be precise, the two models considered in [35] can be written as follows:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =b{r}_S(S,I)S\left(1-\frac{S+I}{\theta }\right)-\beta SW-hS+gI,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =b{r}_I(S,I)I\left(1-\frac{S+I}{\theta }\right)+\beta SW-(\alpha +h+g)I,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \frac{dU}{dt}& =a(U+W)\left(1-\frac{U+W}{K(S+I)}\right)-\gamma IU-(\mu +c)U,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \frac{dW}{dt}& =\gamma IU-(\mu +c)W,\hfill \end{array}$$

(10)

with $r_S(S,I)$ and $r_I(S,I)$ considered to unify the models, referred to as the “frequency-replanting model” and “abundance-replanting model”, and defined as

$$\begin{array}{c}\hfill r_S(S,I)=\left\{\begin{array}{cc}1,\hfill & \text{abundance-replanting }\mathrm{model},\hfill \\ {\displaystyle \frac{S}{S+ϵI},}\hfill & \text{frequency-replanting }\mathrm{model},\hfill \end{array}\right.\end{array}$$

(11)

$$\begin{array}{c}\hfill r_I(S,I)=\left\{\begin{array}{cc}\widehat{ϵ},\hfill & \text{abundance-replanting }\mathrm{model},\hfill \\ {\displaystyle \frac{ϵI}{S+ϵI},}\hfill & \text{frequency-replanting }\mathrm{model}.\hfill \end{array}\right.\end{array}$$

(12)

The parameters $b,\theta ,\beta ,h,g,\alpha ,a,\gamma ,\mu ,c$, and $ϵ$ are positive constants denoting some rates as described in

Table 1. Parameter definitions for system (7)–(12).

Parameter	Definition	Dimension
Host
b	Maximum replanting rate	per day
$\theta$	Maximum plant population pressure	plants per m${}^2$
$\beta$	Vector to host inoculation rate	per vector per day
h	Harvesting rate	per day
g	Recovery rate from infection	per day
$\alpha$	Roguing rate (remotion of infected plants)	per day
$ϵ$	Parameter in frequency-replanting model	-
$\widehat{ϵ}$	Parameter in abundance-replanting model	-
Vector
a	Maximum vector birth rate	per day
$\mu$	Natural vector mortality	per day
c	Insecticide vector mortality	per day
K	Maximum host abundance	vectors per m${}^2$
Pathogen
$\gamma$	Host to vector acquisition rate	per plant per day

. We remark that the parameters $ϵ$/$\widehat{ϵ}$ used in both models have different meanings: In the case of the frequency-replanting model, the selection is based on the composition of susceptible and infected plants, and in the case of the abundance-replanting model, the selection is based on the abundance of infected plants. Moreover, some important remarks about the assumptions used to deduce the system (7)–(12) are the following: (i) plants are propagated vegetatively, by continuous replanting, using plant cuttings, within the crop field; (ii) it is not possible to identify infected plants prior to replanting, so “frequency-of-selection” or “replanting-by-abundance” replanting includes healthy plants and some infected plants; (iii) the growth of both populations is limited by the carrying capacity of the crop field; (iv) the rate of replanting, in any form, is density dependent; (v) there is no latency period in the plant population; and (vi) vectors are infectious throughout their lives, so there is no latency period.

In this paper, we are interested in the analysis of periodic solutions of the plant virus system (7)–(12). We then propose a generalization by assuming that some parameters are positive, continuous, and periodic functions. To be more precise, we present the mathematical model as follows:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =b\left(t\right)r_S(S,I)S\left(1-\frac{S+I}{\theta }\right)-\beta \left(t\right)SW-h\left(t\right)S+g\left(t\right)I,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =b\left(t\right)r_I(S,I)I\left(1-\frac{S+I}{\theta }\right)+\beta \left(t\right)SW-(\alpha +h+g)\left(t\right)I,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \frac{dU}{dt}& =a\left(t\right)(U+W)\left(1-\frac{U+W}{K(S+I)}\right)-\gamma \left(t\right)IU-(\mu +c)\left(t\right)U,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \frac{dW}{dt}& =\gamma \left(t\right)IU-(\mu +c)\left(t\right)W.\hfill \end{array}$$

(16)

Moreover, we consider that the system is supplemented with the following initial condition:

$$\begin{array}{c}\hfill (S,I,U,W)\left(0\right)=(S_0,I_0,U_0,W_0).\end{array}$$

(17)

Figure 1 shows the relationship between the compartments of the model equations.

Theorem 1. 

Assume that $\widehat{ϵ},ϵ,\theta$, and K are positive constants with $\widehat{ϵ}<1$; the initial condition $(S_0,I_0,U_0,W_0)\in {\mathbb{R}}_{+}^4$; and the functions $b$, $\beta$, $h$, $g$, $\alpha$, $\gamma$, μ, and c are positive, continuous, ω-periodic on $[0,\omega ]$. Then, system (13)–(17) has at least one positive ω-periodic solution. Study of a Preliminary Model shown in Appendix A.

A previous result related to Theorem 1 for frequency-replanting models and abundance-replanting models when the coefficients are constant functions was deduced in [35], where the authors achieved the result of existence for non-negative solutions by considering the assumption that $h+\alpha \le \mu +c$ (see [35], theorem A.1). Thus, Theorem 1 is an improvement on the result from [35] since we do not consider the restrictions on $h+\alpha$ and $\mu +c$.

In this paper, we apply topological degree theory, and the organization is as follows. In Section 2, for a self-contained presentation, we recall some notation, concepts, and results related to the statement of Mawhin’s continuation theorem [36]. Moreover, in Section 2 we recall a useful result on the lower estimates for a linear differential equation. In Section 3, we define some particular Banach spaces $X$, Y and operators L and N, and prove that these particular operators satisfy the assumptions of Mawhin’s continuation theorem. In Section 4, we develop the proof of Theorem 1. Finally, in Section 5 and Section 6 we present some numerical examples and the conclusions, respectively.

2. Preliminaries

We present some terminology notation, terminology of topological degree theory, and results that are relevant to proving Theorem 1 by application of Mawhin’s continuation theorem [36].

2.1. Topological Degree Theory Notation an Results

Definition 1. 

Let $L:\mathrm{Dom}L\subset X\to Y$ be a linear operator defined on the normed vector spaces X and Y. If $\mathrm{Im}L$ and $\mathrm{Ker}L$ are such that

$$\begin{array}{c}\hfill \dim \left(\mathrm{Ker}L\right)=\mathrm{codim}\left(\mathrm{Im}L\right)<\infty ,\mathrm{Im}L\mathit{is}\mathit{closed}\mathit{in}Y,\end{array}$$

(18)

are satisfied, L is a Fredholm operator of index zero.

Proposition 1. 

Let L be a Fredholm mapping of index zero from X to Y, then

(i) 

There exist $P:X\to X$ and $Q:Y\to Y$ such that P and Q are continuous projectors satisfying $\mathrm{Im}P=\mathrm{Ker}L$ and $\mathrm{Im}L=\mathrm{Ker}Q=\mathrm{Im}(I-Q)$.

(ii) 

The restriction of L to operator $\mathrm{Dom}L\cap \mathrm{Ker}P$ defines an invertible operator $L_P$, i.e., $L_P{:=L|}_{\mathrm{Dom}L\cap \mathrm{Ker}P}:(I-P)X\to \mathrm{Im}L$ and $K_P=L_P^{-1}$ is well defined.

(iii) 

The sets $\mathrm{Im}Q$ and $\mathrm{Ker}L$ are isomorphic, there exists $J:\mathrm{Im}Q\to \mathrm{Ker}L$ such that J is an isomorphism.

Definition 2. 

Let X and Y be two normed vector spaces and $\Omega \subset X$ be an open bounded set; and L a Fredholm mapping of index zero from X to Y with $P,Q,$ and $K_P$ the operators satisfying the assertions (i)–(iii) of Proposition 1. The operator $N:X\to Y$ is called L-compact on $\overline{\Omega }$ if the following assertions are satisfied: (a) N is continuous, (b) $QN\left(\overline{\Omega }\right)$ is bounded, and (c) the operator $K_P(I-Q)N$ is compact on the set $\overline{\Omega }$.

Definition 3. 

Let us consider $\Omega \subset {\mathbb{R}}^n,f:\Omega \to {\mathbb{R}}^n,y\in {\mathbb{R}}^n$ and $N_f=\{x\in \Omega :J_f\left(x\right)=0\}$, with $J_f\left(x\right)$ the Jacobian of f at x, such that Ω is open and bounded, $f\in C^1(\Omega ,{\mathbb{R}}^n)\cap C(\overline{\Omega },{\mathbb{R}}^n)$, and $y\in {\mathbb{R}}^n\setminus f(\partial \Omega \cup N_f)$. Then, $\mathrm{deg}\{f,\Omega ,y\}={\sum }_{x\in f^{-1}\left(y\right)}\mathrm{sgn}{J}_f\left(x\right),$ defines the degree $\mathrm{deg}\{f,\Omega ,y\}$, with the agreement ${\displaystyle \sum _{\varnothing }=0}$.

Theorem 2. 

Assume that X and Y are Banach spaces, $\Omega \subset X$ is open bounded, L is a Fredholm mapping of index zero from X to Y, and N is an L-compact operator on $\overline{\Omega }$. If the following hypotheses are met:

(C${}_1$) 

For each $(\lambda ,x)\in (0,1)\times (\partial \Omega \cap \mathrm{Dom}L)$ the relation $Lx\ne \lambda Nx$ is satisfied; for each $(\lambda ,x)\in (0,1)\times (\partial \Omega \cap \mathrm{Dom}L)$.

(C${}_2$) 

For each $x\in \partial \Omega \cap \mathrm{Ker}L$ the relation $QNx\ne 0$ is satisfied;

(C${}_3$) 

$\mathrm{deg}(JQN,\Omega \cap \mathrm{Ker}L,0)\ne 0$.

Then, there exists at least one $x\in \mathrm{Dom}L\cap \overline{\Omega }$ satisfying the operator equation $Lx=Nx$.

2.2. A Useful Auxiliary Result

Proposition 2 

([37]). Let $\psi :[0,\omega ]\subset {\mathbb{R}}^{+}\to \mathbb{R}$ be a function such that it is absolutely continuous and satisfies the differential inequality

$$\begin{array}{c}\hfill \frac{d}{dt}\psi \left(t\right)+m\left(t\right)\psi \left(t\right)\ge 0,\forall t\in [0,\omega ],\end{array}$$

(19)

for a function $m\in L^1\left([0,\omega ]\right)$ which is bounded on $[0,\omega ]$, i.e., there exists $m_1,m_2\in {\mathbb{R}}_{+}$ such that $0<m_1\le m\left(t\right)\le m_2$ for all $t\in [0,\omega ]$. Then, the assumption $\psi \left(0\right)>0$ implies the estimate $\psi \left(t\right)\ge \psi \left(0\right)exp(-m_2\omega )>0$ for all $t\in [0,\omega ].$

3. The Banach Space of Continuous Periodic Functions and Particular Operators L and N

3.1. Definitions of Spaces X and Y and the Operators L and N

We consider that the spaces X and Y are coincident and defined as the Banach space of continuous and $\omega$-periodic functions. More precisely:

$$\begin{array}{cc}\hfill & X=Y=\{\mathbf{x}={(x_1,x_2,x_3,x_4)}^T\in C(\mathbb{R};{\mathbb{R}}^4):\hfill \\ \hfill & x_i(t+\omega )=x_i\left(t\right),i=\overline{1,4}∥\mathbf{x}∥:=\sum _{i=1}^4\underset{t\in [0,\omega ]}{max}\left|x_i\left(t\right)\right|<\infty \}.\hfill \end{array}$$

(20)

The operators L and N from X to Y are defined as follows:

$$\begin{array}{cc}\hfill & L\left({\mathbf{x}}^T\right)={\left(\frac{d\mathbf{x}}{dt}\right)}^T,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & N\left({\mathbf{x}}^T\right)=(b\left(t\right)r_S(exp\left(x_1\right),exp\left(x_2\right))\left[1-\frac{exp\left(x_1\right)+exp\left(x_2\right)}{\theta }\right]-\beta \left(t\right)exp\left(x_4\right)-h\left(t\right)\hfill \\ \hfill & +g\left(t\right)exp(x_2-x_1),b\left(t\right)r_I(exp\left(x_1\right),exp\left(x_2\right))\left[1-\frac{exp\left(x_1\right)+exp\left(x_2\right)}{\theta }\right]\hfill \\ \hfill & +\beta \left(t\right)exp(x_1+x_4-x_2)-(\alpha +h+g)\left(t\right),\hfill \\ \hfill & a\left(t\right)(1+exp(x_4-x_3))\left[1-\frac{exp\left(x_3\right)+exp\left(x_4\right)}{K(exp\left(x_1\right)+exp\left(x_2\right))}\right]-\gamma \left(t\right)exp\left(x_2\right)-(\mu +c)\left(t\right),\hfill \\ \hfill & \gamma \left(t\right)exp(x_2+x_3-x_4)-(\mu +c)\left(t\right){)}^T.\hfill \end{array}$$

(22)

Moreover, we consider $\Omega \subset X$ the open ball of X centered at $\mathbf{0}\in X$ and radius $\sigma$, i.e., $\Omega$ is the following set:

$$\begin{array}{c}\hfill \Omega =\left\{\mathbf{x}\in X:\parallel \mathbf{x}\parallel <\sigma ,\sigma \in {\mathbb{R}}^{+}\right\}.\end{array}$$

(23)

Particularly, in some of the results given below, we fix $\sigma$ by defining it in terms of the coefficients.

3.2. L Is a Fredholm Operator of Index Zero Operator and N Is an L-Compact Operator

Lemma 1. 

Let us consider the operators L and N give by the relations on (21) and (22) with $X=Y$ defined in (20) and Ω given by the set in (23). Then, the following assertions are satisfied:

(a) 

L is a Fredholm operator of index zero and

$$\begin{array}{c}\hfill \mathrm{Ker}L\cong {\mathbb{R}}^4,\mathrm{Im}L=\left\{\mathbf{y}\in Y:{\int }_0^{\omega }\mathbf{y}{\left(\tau \right)}^{T}d\tau =0\right\}.\end{array}$$

(24)

(b) 

N is a continuous operator.

(c) 

The operators P and Q defined by

$$\begin{array}{c}\hfill P:X\to X,Q:Y\to Y,P\left({\mathbf{x}}^T\right)=Q\left({\mathbf{x}}^T\right)=\frac{1}{\omega }{\int }_0^{\omega }\mathbf{x}{\left(\tau \right)}^{T}d\tau \end{array}$$

(25)

satisfy the relations (i)–(iii) given in Proposition 1.

(d) 

If the assumptions of Theorem 1 are met, the operator N is L-compact on $\overline{\Omega }$.

Proof. 

Proof of (a). Utilizing the definition of L, we deduce the following:

$$\begin{array}{cc}\hfill \mathrm{Ker}\left(L\right)& =\left\{{\mathbf{x}}^T\in \mathrm{Dom}\left(L\right):L\left(\mathbf{x}\right)={(0,0,0,0)}^T\right\}\hfill \\ \hfill & =\left\{{\mathbf{x}}^T\in \mathrm{Dom}\left(L\right):\frac{d{x}_i}{dt}\left(t\right)=0,i=\overline{1,4}\right\}\hfill \\ \hfill & =\left\{{\mathbf{x}}^T\in \mathrm{Dom}\left(L\right):\mathbf{x}\left(t\right)=\mathbf{x}\left(0\right),t\in [0,\omega ]\right\}\cong {\mathbb{R}}^4.\hfill \end{array}$$

Consequently, $\dim \left(\mathrm{Ker}\right(L\left)\right)=4<\infty$. Moreover, employing a result from elementary algebra, it is established that there exists a unique isomorphism $\phi :X/\mathrm{Ker}\left(L\right)\to \mathrm{Im}\left(L\right),$ which ensures the commutativity visualized by the following scheme:

Here, $\pi :X\to X/\mathrm{Ker}\left(L\right)$ represents the canonical projection, and $i:\mathrm{Im}\left(L\right)\to Y$ signifies the inclusion. This implies that $X/\mathrm{Ker}\left(L\right)\cong \mathrm{Im}\left(L\right).$ Additionally, based on elementary algebra, we can identify the isomorphisms $X\cong \mathrm{Ker}\left(L\right)⨁X/\mathrm{Ker}\left(L\right)$ and $Y\cong \mathrm{Im}\left(L\right)⨁X/\mathrm{Im}\left(L\right)$. It follows that $X\cong Ker\left(L\right)⨁X/ker\left(L\right)\cong \mathrm{Ker}\left(L\right)⨁\mathrm{Im}\left(L\right)$ Given that $\mathrm{Im}\left(L\right)⨁X/\mathrm{Im}\left(L\right)\cong Y=X\cong \mathrm{Ker}\left(L\right)⨁\mathrm{Im}\left(L\right)$, this implies that $\mathrm{Ker}\left(L\right)\cong X/\mathrm{Im}\left(L\right)$, and thus, $\dim \left(\mathrm{Ker}\right(L\left)\right)=\mathrm{cod}\left(\mathrm{Im}\right(L\left)\right)=4<\infty$. Therefore, by Definition 1 we have that L is a Fredholm operator of index zero.

We have deduced that $\mathrm{Ker}L\cong {\mathbb{R}}^4$. To characterize the set defining $\mathrm{Im}\left(L\right)$, we choose arbitrarily ${\mathbf{y}}^T\in \mathrm{Im}\left(L\right)$, and find $\mathbf{x}\in \mathrm{Dom}\left(L\right)$ such that $L{\mathbf{x}}^T={\mathbf{y}}^T$, i.e., ${\mathbf{x}}^{\prime }\left(t\right)=\mathbf{y}\left(t\right)$. The periodicity of $\mathbf{x}$, implies that ${\int }_0^{\omega }{\mathbf{x}}^{\prime }\left(\tau \right)d\tau =0$. Consequently, ${\mathbf{y}}^T\in \mathrm{Im}\left(L\right)$ if and only if ${\int }_0^{\omega }\mathbf{y}\left(\tau \right)d\tau =0$, and thus, $\mathrm{Im}\left(L\right)$ is characterized by the set defined in (24).

Proof of (b). Consider the sequence $\left\{{\mathbf{x}}_n\right\}\subset X$ and $\overline{\mathbf{x}}\in X$ such that ${\mathbf{x}}_n\to \overline{\mathbf{x}}$ in the topology induced by the norm of X. Based on the definition of the operator N as given in (22), and applying the inequality component-wise,

$$\begin{array}{c}\hfill |exp\left(z_2\right)-exp\left(z_1\right)|=\left|{\int }_{z_2}^{z_1}exp\left(s\right)ds\right|\le max\left\{exp\left(z_1\right),exp\left(z_2\right)\right\}|z_2-z_1|,\forall z_1,z_2\in \mathbb{R},\end{array}$$

we establish the existence of a constant $C>0$ depending only on $h$, $g$, $\alpha$, $\gamma$, $\mu$, $c$, $ϵ$, $\theta$, and K, such that $\parallel N\left({\mathbf{x}}_n\right)-N\left(\overline{\mathbf{x}}\right)\parallel \le C\parallel {\mathbf{x}}_n-\overline{\mathbf{x}}\parallel .$ This leads to the conclusion that $N\left({\mathbf{x}}_n\right)\to N\left(\overline{\mathbf{x}}\right)$ in the norm sense, proving that N is continuous.

Proof of (c). Our proof of $\mathrm{Ker}L=\mathrm{Im}P$ is developed by double inclusion. Firstly, we obtain that $\mathrm{Ker}L\subset \mathrm{Im}P$ as follows: from the isomorphism $\mathrm{Ker}L\cong {\mathbb{R}}^4$ (see Lemma 1-(a)), we note that ${\mathbf{x}}^T\in \mathrm{Ker}L$ is equivalent to the fact that $\mathbf{x}\left(t\right)$ is a constant function on $[t_0,\infty ]$, which implies that $\mathbf{x}\in \mathrm{Im}P$, since $P\left({\mathbf{x}}^T\right)={\mathbf{x}}^T$ for constant $\mathbf{x}\left(t\right)$. Secondly, we show that $\mathrm{Im}P\subset \mathrm{Ker}L$: for ${\mathbf{y}}^T\in \mathrm{Im}P$, there exists $\mathbf{z}\in X$ such that $P\left({\mathbf{z}}^T\right)={\mathbf{y}}^T$, and from Equation (25), we have ${\omega }^{-1}{\int }_0^{\omega }\mathbf{z}{\left(\tau \right)}^{T}d\tau ={\mathbf{y}}^T$, which implies, by differentiation, that $L\left({\mathbf{y}}^T\right)=0$ or $\mathbf{y}\in \mathrm{Ker}L.$

We prove that $\mathrm{Ker}Q=\mathrm{Im}L$ by using the definition of the operator Q provided in (25). We observe the fact that ${\mathbf{y}}^T\in \mathrm{Ker}Q$ is rewritten as ${\int }_0^{\omega }\mathbf{y}{\left(\tau \right)}^{T}d\tau =0$, which, according to (24), is equivalent to ${\mathbf{y}}^T\in \mathrm{Im}L.$

The proof that $\mathrm{Im}(I-Q)=\mathrm{Im}L$ is carried out through double inclusion. Given ${\mathbf{y}}^T\in \mathrm{Im}(I-Q)$, there exists $\mathbf{z}\in X$ such that $(I-Q)\left({\mathbf{z}}^T\right)={\mathbf{y}}^T$, implying

$$\begin{array}{c}\hfill {\int }_0^{\omega }\mathbf{y}{\left(\tau \right)}^{T}d\tau ={\int }_0^{\omega }\left(\mathbf{z}{\left(\tau \right)}^T-\frac{1}{\omega }{\int }_0^{\omega }\mathbf{z}{\left(m\right)}^{T}dm\right)d\tau =(0,0,0,0),\end{array}$$

and, according to (24), we obtain that $\mathbf{y}{\left(\tau \right)}^T\in \mathrm{Im}L.$ Hence, we deduce that $\mathrm{Im}(I-Q)\subset \mathrm{Im}L$. Similarly, we establish the inclusion $\mathrm{Im}L\subset \mathrm{Im}(I-Q)$.

Let us prove the existence of the operators $K_P$ and $L_P.$ The operator $L_P$ is the restriction of L to $\mathrm{Dom}L\cap \mathrm{Ker}P,$ i.e., $L_P:\mathrm{Dom}L\cap \mathrm{Ker}P\to \mathrm{Im}L$ and $L_P=L$ on $\mathrm{Dom}L\cap \mathrm{Ker}P$. The operator $K_P$ is the inverse operator of $L_P$, defined by the relation

$$\begin{array}{cc}\hfill K_P& \left({\mathbf{x}}^T\right)\left(t\right)={\int }_0^t\mathbf{x}{\left(\tau \right)}^{T}d\tau -\frac{1}{\omega }{\int }_0^{\omega }{\int }_0^{\eta }\mathbf{x}{\left(m\right)}^{T}dmd\eta .\hfill \end{array}$$

(26)

Then, to prove that $K_P=L_P^{-1}$, we apply the following identity:

$$\begin{array}{ccc}& & {\int }_0^t\frac{d}{ds}\mathbf{x}\left(s\right)ds-\frac{1}{\omega }{\int }_0^{\omega }{\int }_0^t\frac{d}{dm}\mathbf{x}\left(m\right)dmdt=\mathbf{x}\left(t\right),\hfill \end{array}$$

which is satisfied only on $\mathrm{Dom}L\cap \mathrm{Ker}P.$

Proof of (d). The proof is facilitated by applying Definition 2. Considering the definition of the operator Q, we note that $QN\left({\mathbf{x}}^T\right)=\frac{1}{\omega }{\int }_0^{\omega }N{\left(\tau \right)}^{T}d\tau .$ Then, for $\mathbf{x}\in \overline{\Omega }$ we obtain that $\parallel QN\left({\mathbf{x}}^T\right)\parallel \le \frac{1}{\omega }{\int }_0^{\omega }\parallel N\parallel d\tau =\parallel N\parallel ,$ which implies that $QN\left(\overline{\Omega }\right)$ is bounded on $\overline{\Omega }$. From the definitions of the operators $K_P,N$, and Q, we have

$$(K_P(I-Q)N)\left({\mathbf{x}}^T\right)\left(t\right)={\int }_0^{t}N{\left(\tau \right)}^{T}d\tau +\left(\frac{1}{2}-\frac{t}{\omega }\right){\int }_0^{\omega }N{\left(\tau \right)}^{T}d\tau -\frac{1}{\omega }{\int }_0^{\omega }{\int }_0^{\eta }N{\left(m\right)}^{T}dmd\eta .$$

Consequently, we deduce the estimate $\parallel K_P(I-Q)N\parallel \le 2\omega \parallel N\parallel ,$ indicating that $\left(K_P(I-Q)N\right)\left(\overline{\Omega }\right)$ is bounded on $\overline{\Omega }$, as N is bounded on $\overline{\Omega }$. Additionally, we clearly deduce that

$$|\left(K_P(I-Q)N\right)\left({\mathbf{x}}^T\right)\left(t\right)-\left(K_P(I-Q)N\right)\left({\mathbf{x}}^T\right)\left(s\right)|\le 2\parallel N\parallel |t-s|,\forall t,s\in [t_0,\infty ],$$

or equivalently, $K_P(I-Q)N$ is an equicontinuous operator. By applying the Arzela–Ascoli theorem, it follows that $K_P(I-Q)N$ is a compact operator on $\overline{\Omega }$. Hence, we have satisfied the conditions of Definition 2, thereby concluding that N is L-compact on $\overline{\Omega }$. □

3.3. Results Related to the Operator Equation $L\mathbf{x}=\lambda N\mathbf{x}$ for $(\lambda ,\mathbf{x})\in (0,1)\times (\partial \Omega \cap DomL)$

Before introducing the results, we define some notation. For any function $f:[0,\omega ]\to \mathbb{R}$ continuous, positive, and bounded on $[0,\omega ]$, we consider

$$\begin{array}{c}\hfill \overline{f}=\frac{1}{\omega }{\int }_0^{\omega }f\left(t\right)dt,f^{\perp }=\underset{t\in [0,\omega ]}{min}f\left(t\right),f^{\top }=\underset{t\in [0,\omega ]}{max}f\left(t\right).\end{array}$$

(27)

Proposition 3. 

Assuming the conditions of Lemma 1 and Theorem 1 are satisfied, and considering

$$\begin{array}{cc}\hfill & {\mathsf{\Psi }}_1=2exp(\sigma +\omega b^{\top }),{\mathsf{\Psi }}_2={\mathsf{\Psi }}_1,{\mathsf{\Psi }}_3=2exp(\sigma +\omega a^{\top }),{\mathsf{\Psi }}_3={\mathsf{\Psi }}_4,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & {\Lambda }_1=exp\left(x_1\left(0\right)\right)exp\left(-\omega \left[2{\theta }^{-1}{b}^{\top }{\mathsf{\Psi }}_1+{\beta }^{\top }{\mathsf{\Psi }}_4+h^{\top }\right]\right),\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & {\Lambda }_2=exp\left(x_2\left(0\right)\right)exp\left(-\omega \left[2{\theta }^{-1}{b}^{\top }{\mathsf{\Psi }}_1+{(\alpha +h+g)}^{\top }\right]\right),\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & {\Lambda }_3=exp\left(x_4\left(0\right)\right)exp\left(-\omega \left[2k^{-1}{a}^{\top }{\mathsf{\Psi }}_1{({\Lambda }_1+{\Lambda }_2)}^{-1}+{\gamma }^{\top }+{(\mu +c)}^{\top }\right]\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & {\Lambda }_4=exp\left(x_4\left(0\right)\right)exp\left(-\omega {(\mu +c)}^{\top }\right).\hfill \end{array}$$

(32)

If $(\lambda ,\mathbf{x})\in (0,1)\times (\partial \Omega \cap \mathrm{Dom}L)$ is a solution of the operator equation $L\mathbf{x}=\lambda N\mathbf{x}$, then the following estimates

$$\begin{array}{c}\hfill 0<{\Lambda }_i\le exp\left(x_i\left(t\right)\right)\le {\mathsf{\Psi }}_i,i=\overline{1,4},\end{array}$$

(33)

are satisfied for $t\in [0,\omega ]$.

Proof. 

Given the definitions of the operators L and N, we deduce that the operator equation $L\mathbf{x}=\lambda N\mathbf{x}$ is equivalent to the following ordinary differential equations system:

$$\begin{array}{cc}\hfill & \frac{d{x}_1}{dt}=\lambda [b\left(t\right)r_S\left(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right)\right)\left(1-\frac{exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)}{\theta }\right)\hfill \\ \hfill & -\beta \left(t\right)exp\left(x_4\left(t\right)\right)-h\left(t\right)+g\left(t\right)exp(x_2\left(t\right)-x_1\left(t\right))],\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & \frac{d{x}_2}{dt}=\lambda [b\left(t\right)r_I(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right))\left(1-\frac{exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)}{\theta }\right)\hfill \\ \hfill & +\beta \left(t\right)exp(x_1\left(t\right)+x_4\left(t\right)-x_2\left(t\right))-(\alpha +h+g)\left(t\right)],\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill & \frac{d{x}_3}{dt}=\lambda [a\left(t\right)(1+exp(x_4\left(t\right)-x_3\left(t\right)))\left(1-\frac{exp\left(x_3\left(t\right)\right)+exp\left(x_4\left(t\right)\right)}{K(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right))}\right)\hfill \\ \hfill & -\gamma \left(t\right)exp\left(x_2\left(t\right)\right)-(\mu +c)\left(t\right)],\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & \frac{d{x}_4}{dt}=\lambda \left[\gamma \left(t\right)exp(x_2\left(t\right)+x_3\left(t\right)-x_4\left(t\right))-(\mu +c)\left(t\right)\right].\hfill \end{array}$$

(37)

By multiplying (34) by $exp\left(x_1\left(t\right)\right)$, (35) by $exp\left(x_2\left(t\right)\right)$, (36) by $exp\left(x_3\left(t\right)\right)$, (37) by $exp\left(x_3\left(t\right)\right)$, and rearranging the terms, we deduce that

$$\begin{array}{cc}\hfill & \frac{dexp\left(x_1\right)}{dt}+\frac{\lambda b\left(t\right)}{\theta }{r}_S\left(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right)\right)exp\left(x_1\left(t\right)\right)\left(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)\right)\hfill \\ \hfill & +\lambda \beta \left(t\right)exp(x_4\left(t\right)+x_1\left(t\right))+h\left(t\right)exp\left(x_1\left(t\right)\right)\hfill \\ \hfill & =\lambda b\left(t\right)r_S\left(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right)\right)exp\left(x_1\left(t\right)\right)+\lambda g\left(t\right)exp\left(x_2\left(t\right)\right),\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill & \frac{dexp\left(x_2\right)}{dt}+\frac{\lambda b\left(t\right)}{\theta }{r}_I\left(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right)\right)exp\left(x_2\left(t\right)\right)\left(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)\right)\hfill \\ \hfill & +\lambda (\alpha +h+g)\left(t\right)exp\left(x_2\left(t\right)\right)\hfill \\ \hfill & =\frac{\lambda b\left(t\right)}{\theta }b\left(t\right)r_I(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right))+\lambda \beta \left(t\right)exp(x_1\left(t\right)+x_4\left(t\right)),\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill & \frac{d\exp \left(x_3\right)}{dt}+\lambda a\left(t\right)\frac{{\left(exp\left(x_3\left(t\right)\right)+exp\left(x_4\left(t\right)\right)\right)}^2}{K(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right))}+\lambda \gamma \left(t\right)exp(x_2\left(t\right)+x_3\left(t\right))\hfill \\ \hfill & +\lambda (\mu +c)\left(t\right)exp\left(x_3\left(t\right)\right)\hfill \\ \hfill & =\lambda a\left(t\right)(exp\left(x_3\left(t\right)\right)+exp\left(x_4\left(t\right)\right)),\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill & \frac{dexp\left(x_4\right)}{dt}+\lambda (\mu +c)\left(t\right)exp\left(x_4\left(t\right)\right)=\lambda \gamma \left(t\right)exp(x_2\left(t\right)+x_3\left(t\right)).\hfill \end{array}$$

(41)

In order to obtain the upper bounds in (33), we apply the Gronwall inequality as follows. By adding (38) and (39), rearranging the terms, and using the fact that the definitions of $r_S$ and $r_I$ provided in (12) imply that the relations $r_S(S,I)\in (0,1]$ and $r_I(S,I)\in (0,1]$ are satisfied for any $(S,I)\in {\mathbb{R}}_{+}^2$, we deduce that

$$\begin{array}{cc}\hfill & \frac{d}{dt}\left(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)\right)\hfill \\ \hfill & +\lambda b\left(t\right)\left(r_S(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right))exp\left(x_1\left(t\right)\right)+r_I(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right))exp\left(x_2\left(t\right)\right)\right)\hfill \\ \hfill & \times \left(\frac{exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)}{\theta }\right)+\lambda h\left(t\right)\left(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right))+\lambda \alpha \left(t\right)exp(x_2\left(t\right)\right)\hfill \\ \hfill & =\lambda b\left(t\right)\left(r_S(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right))exp\left(x_1\left(t\right)\right)+r_I(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right))exp\left(x_2\left(t\right)\right)\right)\hfill \\ \hfill & \le \lambda b\left(t\right)\left(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)\right).\hfill \end{array}$$

(42)

Similarly, by adding (40) and (41), and rearranging the terms, we obtain

$$\begin{array}{cc}\hfill & \frac{d}{dt}\left(exp\left(x_3\left(t\right)\right)+exp\left(x_4\left(t\right)\right)\right)+\lambda a\left(t\right)\left(\frac{{\left(exp\left(x_3\left(t\right)\right)+exp\left(x_4\left(t\right)\right)\right)}^2}{K(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right))}\right)\hfill \\ \hfill & +\lambda (\mu +c)\left(t\right)exp\left(x_3\left(t\right)\right)=\lambda a\left(t\right)\left(exp\left(x_3\left(t\right)\right)+exp\left(x_4\left(t\right)\right)\right)\hfill \end{array}$$

(43)

Applying the Gronwall inequality to the relations (42) and (43), we deduce the inequalities

$$\begin{array}{cc}\hfill & exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)\le \left(exp\left(x_1\left(0\right)\right)+exp\left(x_2\left(0\right)\right)\right)exp\left({\int }_0^t\lambda b\left(\tau \right)d\tau \right),\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill & exp\left(x_3\left(t\right)\right)+exp\left(x_4\left(t\right)\right)\le \left(exp\left(x_3\left(0\right)\right)+exp\left(x_4\left(0\right)\right)\right)exp\left({\int }_0^t\lambda a\left(\tau \right)d\tau \right).\hfill \end{array}$$

(45)

Then, using the fact that $\parallel \mathbf{x}\left(0\right)\parallel =\sigma$ and $\lambda \in (0,1)$, we obtain the upper bounds in the estimate (33).

On the other hand, to obtain the lower bounds in the estimate (33), we apply Proposition 2. From (38)–(41), we deduce that the inequalities

$$\begin{array}{cc}\hfill & \frac{dexp\left(x_1\right)}{dt}+M_1\left(t\right)exp\left(x_1\left(t\right)\right)\hfill \\ \hfill & \ge \lambda b\left(t\right)r_S\left(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right)\right)exp\left(x_1\left(t\right)\right)+g\left(t\right)exp\left(x_2\left(t\right)\right)\ge 0,\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill & \frac{dexp\left(x_2\right)}{dt}+M_2\left(t\right)exp\left(x_2\left(t\right)\right)\hfill \\ \hfill & \ge \frac{\lambda b\left(t\right)}{\theta }b\left(t\right)r_I(exp\left(x_1\left(t\right)\right),exp\left(x_2\left(t\right)\right))+\lambda \beta \left(t\right)exp(x_1\left(t\right)+x_4\left(t\right))\ge 0,\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill & \frac{dexp\left(x_3\right)}{dt}+M_3\left(t\right)exp\left(x_3\left(t\right)\right)\ge \lambda a\left(t\right)(exp\left(x_3\left(t\right)\right)+exp\left(x_4\left(t\right)\right))\ge 0,\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill & \frac{dexp\left(x_4\right)}{dt}+M_4\left(t\right)exp\left(x_4\left(t\right)\right)=\lambda \gamma \left(t\right)exp(x_2\left(t\right)+x_3\left(t\right))\ge 0,\hfill \end{array}$$

(49)

are satisfied for $t\in [0,\omega ],$ where the functions $M_i,i=1,\dots ,4$ are given by

$$\begin{array}{cc}\hfill M_1\left(t\right)& =\frac{\lambda b\left(t\right)}{\theta }\underset{t\in [0,\omega ]}{max}\left(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)\right)+\lambda \beta \left(t\right)\underset{t\in [0,\omega ]}{max}exp\left(x_4\left(t\right)\right)+\lambda h\left(t\right)\hfill \\ \hfill M_2\left(t\right)& =\frac{\lambda b\left(t\right)}{\theta }\underset{t\in [0,\omega ]}{max}\left(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)\right)+\lambda (\alpha \left(t\right)+h\left(t\right)+g\left(t\right))\hfill \\ \hfill M_3\left(t\right)& =\frac{\lambda a\left(t\right)}{K}\frac{{\displaystyle \underset{t\in [0,\omega ]}{max}\left(exp\left(x_3\left(t\right)\right)+exp\left(x_4\left(t\right)\right)\right)}}{{\displaystyle \underset{t\in [0,\omega ]}{min}\left(exp\left(x_1\left(t\right)\right)+exp\left(x_2\left(t\right)\right)\right)}}+\lambda \gamma \left(t\right)\underset{t\in [0,\omega ]}{max}exp\left(x_2\left(t\right)\right)+M_4\left(t\right),\hfill \\ \hfill M_4\left(t\right)& =\lambda \left(\mu \right(t)+c(t\left)\right).\hfill \end{array}$$

We observe that the upper bounds in (33) imply that

$$\begin{array}{cc}\hfill & \lambda h^{\perp }\le M_1\left(t\right)\le 2{\theta }^{-1}{b}^{\top }{\mathsf{\Psi }}_1+{\beta }^{\top }{\mathsf{\Psi }}_4+h^{\top },t\in [0,\omega ],\hfill \\ \hfill & \lambda {(\alpha +h+g)}^{\perp }\le M_2\left(t\right)\le 2{\theta }^{-1}{b}^{\top }{\mathsf{\Psi }}_1+{(\alpha +h+g)}^{\top },t\in [0,\omega ],\hfill \\ \hfill & \lambda {(\mu +c)}^{\perp }\le M_4\left(t\right)\le {(\mu +c)}^{\top },t\in [0,\omega ].\hfill \end{array}$$

Then, from application of Proposition 2 to (46), (47), and (49), we deduce the lower bounds in (33) for $i=1,2,4$, respectively. Moreover, we notice that

$$\begin{array}{cc}\hfill & \lambda {(\mu +c)}^{\perp }\le M_3\left(t\right)\le 2k^{-1}{a}^{\top }{\mathsf{\Psi }}_3{({\Lambda }_1+{\Lambda }_2)}^{-1}+{\gamma }^{\top }+{(\mu +c)}^{\top },t\in [0,\omega ],\hfill \end{array}$$

and by application of Proposition 2 to (48), we deduce the lower bound given in (33) for $i=3.$ □

Lemma 2. 

Let us consider that the assumptions of Lemma 1 and Theorem 1 are satisfied. Then, there exist $\mathbf{\Psi },\mathbf{\zeta },\mathit{\Lambda }\in {\mathbb{R}}_{+}^4$ such that the estimates

$$\begin{array}{cc}\hfill & x_i\left(t\right)<ln\left({\mathsf{\Psi }}_i\right)+{\zeta }_i,i=\overline{1,4},t\in [0,\omega ],\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill & ln\left({\Lambda }_i\right)<x_i\left(t\right),i=\overline{1,4},t\in [0,\omega ],\hfill \end{array}$$

(51)

are satisfied for any $(\lambda ,\mathbf{x})\in (0,1)\times (\partial \Omega \cap \mathrm{Dom}L)$ that is a solution of the operator equation $L\mathbf{x}=\lambda N\mathbf{x}$.

Proof. 

Let us consider a pair $(\lambda ,\mathbf{x})\in (0,1)\times (\partial \Omega \cap \mathrm{Dom}L)$ that is a solution of the operator equation $L\mathbf{x}=\lambda N\mathbf{x}$. From Proposition (3), we deduce that

$$\begin{array}{c}\hfill 0<\omega {\Lambda }_i={\int }_0^{\omega }{\Lambda }_{i}d\tau \le {\int }_0^{\omega }exp\left(x_i\left(\tau \right)\right)d\tau \le {\int }_0^{\omega }{\mathsf{\Psi }}_{i}d\tau =\omega {\mathsf{\Psi }}_i,i=\overline{1,4}.\end{array}$$

Then, applying the intermediate value theorem for integrals, there exists ${\tau }_i\in [0,\omega ]$ such that $\omega exp\left(x_i\left({\tau }_i\right)\right)<\omega {\mathsf{\Psi }}_i$ for $i=\overline{1,4},$ or equivalently

$$\begin{array}{c}\hfill x_i\left({\tau }_i\right)<ln\left({\mathsf{\Psi }}_i\right),i=\overline{1,4}.\end{array}$$

(52)

Moreover, recalling that the operator equation $L\mathbf{x}=\lambda N\mathbf{x}$ is equivalent to system (34)–(37), we can take the modulus in each equation and integrate it over $[0,\omega ]$. By utilizing the equations in system (34)–(37), and considering the fact that $\lambda \in (0,1)$, we obtain:

$$\begin{array}{c}\hfill {\int }_0^{\omega }\left|\frac{d{x}_i\left(t\right)}{dt}\right|dt\le {\zeta }_i,i=\overline{1,4},\end{array}$$

(53)

where

$$\begin{array}{cc}\hfill {\zeta }_1& =\omega \left[b^{\top }(1+2{\mathsf{\Psi }}_1{\theta }^{-1})+2{\beta }^{\top }{\mathsf{\Psi }}_3+h^{\top }+g^{\top }{\mathsf{\Psi }}_1{\Lambda }_1^{-1}\right],\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill {\zeta }_2& =\omega \left[b^{\top }(1+2{\mathsf{\Psi }}_1{\theta }^{-1})+2{\beta }^{\top }{\mathsf{\Psi }}_1{\mathsf{\Psi }}_2{\Lambda }_2^{-1}+{(\alpha +h+g)}^{\top }\right],\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill {\zeta }_3& =\omega \left[a^{\top }(1+K{\mathsf{\Psi }}_4{\Lambda }_3^{-1})(1+2{\mathsf{\Psi }}_3){({\Lambda }_1+{\Lambda }_2)}^{-1}+{\gamma }^{\top }{\mathsf{\Psi }}_3+{(\mu +c)}^{\top }\right],\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill {\zeta }_4& =\omega \left[2{\gamma }^{\top }{\mathsf{\Psi }}_1{\mathsf{\Psi }}_3{\Lambda }_4^{-1}+{(\mu +c)}^{\top }\right].\hfill \end{array}$$

(57)

From (52) and (53), and the fundamental theorem of calculus, we deduce that

$$\begin{array}{c}\hfill x_i\left(t\right)=x_i\left({\tau }_i\right)+{\int }_{{\tau }_i}^t\frac{d{x}_i\left(s\right)}{dt}ds<ln\left({\mathsf{\Psi }}_i\right)+{\int }_{{\tau }_i}^t\frac{d{x}_i\left(s\right)}{dt}ds<ln\left({\mathsf{\Psi }}_1\right)+{\zeta }_i.\end{array}$$

Consequently, we conclude that (50) is valid. Furthermore, the lower bound given in Proposition (3), stating $exp\left(x_i\left(t\right)\right)\ge {\Lambda }_i$ for $i=\overline{1,4}$ and $t\in [0,\omega ]$, allows us to deduce (51). Thus, the sets $\mathbf{\Psi },\mathbf{\zeta },$ and $\mathit{\Lambda }$, defined by (28)–(32) and (54)–(57), ensure that both (50) and (51) are satisfied. □

Lemma 3. 

Let us consider the assumptions and notation of Lemma 2. Then, the estimates

$$\begin{array}{cc}\hfill & x_i\left(t\right)<ln\left({\mathsf{\Psi }}_i\right),i=\overline{1,4},t\in [0,\omega ],\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill & ln\left({\Lambda }_i\right)<x_i\left(t\right),i=\overline{1,4},t\in [0,\omega ],\hfill \end{array}$$

(59)

are satisfied for any $\mathbf{x}\in Ker\left(L\right)$ that is a solution of the operator equation $QN\left({\mathbf{x}}^T\right)=0.$

Proof. 

From Lemma 1, we know that $\mathrm{Ker}\left(L\right)\cong {\mathbb{R}}^4$, implying that ${\mathbf{x}}^T\in \mathrm{Ker}\left(L\right)$ is equivalent to $\mathbf{x}\left(t\right)=\mathbf{c}$, where $\mathbf{c}$ is a constant vector. Suppose that ${\mathbf{x}}^T$ satisfies the equation $QN\left(\mathbf{x}{\left(t\right)}^T\right)={\mathbf{0}}^T$, or equivalently $QN\left({\mathbf{c}}^T\right)={\mathbf{0}}^T$. From the definition of Q given in (25), it follows that:

$$\begin{array}{cc}\hfill QN\left({\mathbf{c}}^T\right)& ={\mathbf{0}}^T\hfill \\ \hfill \frac{1}{\omega }{\int }_0^{\omega }N\left({\mathbf{c}}^T\right)d\tau & ={\mathbf{0}}^T\hfill \\ \hfill N\left({\mathbf{c}}^T\right)& ={\mathbf{0}}^T\hfill \\ \hfill \lambda N\left({\mathbf{c}}^T\right)& =L\left({\mathbf{c}}^T\right)\hfill \end{array}$$

Thus, for ${\mathbf{x}}^T=\mathbf{c}\in \mathrm{Ker}\left(L\right)$, the operator equations $QN\left({\mathbf{c}}^T\right)={\mathbf{0}}^T$ and $\lambda N\left({\mathbf{c}}^T\right)=L\left({\mathbf{c}}^T\right)$ are equivalent. Moreover, we notice that (53) is satisfied with $\zeta =\mathbf{0}$. Therefore, the proof of (58) and (59) follows as a consequence of Lemma 2. □

3.4. The Items (C1)–(C3) of Theorem 2 Are Satisfied

Lemma 4. 

Let us consider the assumptions and notation of Lemma 2. Additionally, we define the radius σ for the set Ω as follows:

$$\begin{array}{c}\hfill \sigma =max\left\{\sum _{i=1}^4{\left(ln\left({\Lambda }_i\right)\right)}^2,\sum _{i=1}^4{\left(\left|ln\left({\mathsf{\Psi }}_i\right)+{\zeta }_i\right|\right)}^2\right\},\end{array}$$

(60)

where $\mathbf{\Psi },\mathbf{\zeta }$, and $\mathit{\Lambda }$ are as defined in (28)–(32) and (54)–(57). Then, the items (C1)–(C3) of Theorem 2 are satisfied.

Proof. 

The proof of each item is developed basically by applying contradiction arguments.

Proof of (C1). Let us assume that there exists a pair $(\mathbf{x},\lambda )\in (\partial \Omega \cap \mathrm{Dom}(L)\times (0,1)$ such that $L\left({\mathbf{x}}^T\right)=\lambda N\left({\mathbf{x}}^T\right)$. However, according to Lemma 2, this assumption implies that $\mathbf{x}\in Int(\Omega )$, which is a direct contradiction to the initial assumption that $\mathbf{x}\in \partial \Omega .$ Therefore, such a pair $\mathbf{x}\in Int(\Omega )$ cannot exist, and item (C1) is proved.

Proof of (C2). Consider $\mathbf{x}\in \partial \Omega \cap ker\left(L\right)$ such that $QN\left({\mathbf{x}}^T\right)={\mathbf{0}}^T$. From Lemma 3, it follows that $\mathbf{x}\in Int(\Omega )$. This again leads to a contradiction, as it conflicts with our initial assumption that $\mathbf{x}\in \partial \Omega .$

Proof of (C3). We define the application $\Phi :\mathrm{Dom}\left(L\right)\times [0,1]\to X,$ defined as follows:

$$\begin{array}{c}\hfill \Phi ({\mathbf{x}}^T,\epsilon )=\left[\begin{array}{c}\overline{\phi }-\overline{\beta }exp\left(x_4\right)-\overline{h}\hfill \\ \overline{\psi }+\overline{\beta }exp(x_1+x_4-x_2)\hfill \\ \overline{\rho }-\overline{\gamma }exp\left(x_2\right)\hfill \\ \overline{\gamma }exp(x_2+x_3-x_4)-\overline{\eta }\hfill \end{array}\right]+\epsilon \left[\begin{array}{c}\overline{g}exp(x_2-x_1)\hfill \\ -\overline{\alpha +h+g}\hfill \\ \overline{\rho }exp(x_3+x_4)-\overline{\eta }\hfill \\ 0\hfill \end{array}\right].\end{array}$$

Our goal is to show that $\Phi ({\mathbf{x}}^T,\epsilon )\ne 0,$ for $\epsilon \in [0,1]$ and ${\mathbf{x}}^T\in \partial \Omega \cap ker\left(L\right).$ Recall that $ker\left(L\right)\cong {\mathbb{R}}^4,$ such that ${\mathbf{x}}^T$ is a constant function, say ${\mathbf{x}}^T\left(t\right)={(c_1,c_2,c_3,c_4)}^T$. Assuming the existence of $\epsilon \in [0,1]$ and ${(c_1,c_2,c_3,c_4)}^T\in \partial \Omega \cap ker\left(L\right),$ with $∥(c_1,c_2,c_3,c_4)∥=\sigma ,$ such that $\Phi ({(c_1,c_2,c_3,c_4)}^T,\epsilon )=0.$ It follows that

$$\begin{array}{ccc}\hfill 0& =& \overline{\phi }-\overline{\beta }exp\left(c_4\right)-\overline{h}+\epsilon \overline{g}exp(c_2-c_4),\hfill \\ \hfill 0& =& \overline{\psi }+\overline{\beta }exp(c_1+c_4-c_2)-\epsilon \left(\overline{\alpha +h+g}\right),\hfill \\ \hfill 0& =& \overline{\rho }-\overline{\gamma }exp\left(c_2\right)+\epsilon (\overline{\rho }exp(c_4-c_3)-\overline{\eta }),\hfill \\ \hfill 0& =& \overline{\gamma }exp(c_2+c_3)-\overline{\eta }.\hfill \end{array}$$

However, similar arguments to those used in the proof of Lemma 2, imply $∥(c_1,c_2,c_3,c_4)∥<\sigma ,$ contradicting the initial assumption $∥(c_1,c_2,c_3,c_4)∥=\sigma .$

Applying the homotopy invariance theorem of topological degree, performing $J=I:\mathrm{Im}\left(Q\right)\to ker\left(L\right)$ such that ${\mathbf{x}}^T\mapsto {\mathbf{x}}^T,$ using the fact that the system

$$\begin{array}{ccc}\hfill 0& =& \phi -\beta exp\left(x_4\right)-h+gexp(x_2-x_4),\hfill \\ \hfill 0& =& \psi +\beta exp(x_1+x_4-x_2)-(\alpha +h+g),\hfill \\ \hfill 0& =& \rho -\gamma exp\left(x_2\right)+\rho exp(x_4-x_3)-\eta ,\hfill \\ \hfill 0& =& \gamma exp(x_2+x_3)-\eta .\hfill \end{array}$$

has a unique solution ${\mathbf{x}}^{★}\in \partial \Omega \cap ker\left(L\right)$ and by definition 2.3 of [36], we conclude that

$$\begin{array}{cc}\hfill & deg\left(JQN\left({\mathbf{x}}^T\right),\Omega \cap ker\left(L\right),{\mathbf{0}}^T\right)=deg\left(\Phi \left({\mathbf{x}}^T\right),\Omega \cap ker\left(L\right),{\mathbf{0}}^T\right)\hfill \\ \hfill & =deg\left((\overline{\phi }-\overline{\beta }exp\left(x_4\right)-\overline{h},\overline{\psi }+\overline{\beta }exp(x_1+x_4-x_2),\overline{\rho }-\overline{\gamma }exp\left(x_2\right),\overline{\gamma }exp(x_2+x_3-x_4)-\overline{\eta }),\Omega \cap ker\left(L\right),{\mathbf{0}}^T\right)\hfill \\ \hfill & =\mathrm{sgn}\left|\begin{array}{cccc}0& 0& 0& -\overline{\beta }exp\left(x_4\right)\\ \overline{\beta }exp(x_1+x_4-x_2)& -\overline{\beta }exp(x_1+x_4-x_2)& 0& \overline{\beta }exp(x_1+x_4-x_2)\\ 0& -\overline{\gamma }exp\left(x_2\right)& 0& 0\\ 0& \overline{\gamma }exp(x_2+x_3-x_4)& \overline{\gamma }exp(x_2+x_3-x_4)& -\overline{\gamma }exp(x_2+x_3-x_4)\end{array}\right|\hfill \\ \hfill & =\mathrm{sgn}\left(-{\left(\overline{\beta }\overline{\gamma }\right)}^{2}exp(x_1+x_2+x_3+x_4)\right)=-1.\hfill \end{array}$$

Then, $deg\left(JQN\left({\mathbf{x}}^T\right),\Omega \cap ker\left(L\right),{\mathbf{0}}^T\right)\ne 0$. This demonstrates that condition (C3) is valid.

Therefore, for operators L and N, the properties delineated in items (C1), (C2), and (C3) of Theorem 2 are satisfied. □

4. Proof of Theorem 1

4.1. A Change of Variable

Lemma 5. 

Suppose that the ordinary differential equations system (13)–(17) has a solution. Then, $(S,I,U,W)$ satisfies (13)–(17) if and only if $(S^{★},I^{★},U^{★},W^{★}),$ defined as

$$\begin{array}{c}\hfill (S,I,U,W)\left(t\right)=\left(exp\left(S^{★}\left(t\right)\right),exp\left(I^{★}\left(t\right)\right),exp\left(U^{★}\left(t\right)\right),exp\left(W^{★}\left(t\right)\right)\right),\end{array}$$

(61)

is a solution of the following ordinary differential equations system:

$$\begin{array}{cc}\hfill & \frac{d{S}^{★}}{dt}=b_S(exp\left(S^{★}\right),exp\left(I^{★}\right),t)\left(1-\frac{exp\left(S^{★}\right)+exp\left(I^{★}\right)}{\theta }\right)\hfill \\ \hfill & -\beta \left(t\right)exp\left(W^{★}\right)-h\left(t\right)+g\left(t\right)exp(I^{★}-S^{★}),\hfill \end{array}$$

(62a)

$$\begin{array}{cc}\hfill & \frac{d{I}^{★}}{dt}=b_I(exp\left(S^{★}\right),exp\left(I^{★}\right))\left(1-\frac{exp\left(S^{★}\right)+exp\left(I^{★}\right)}{\theta }\right)\hfill \\ \hfill & +\beta \left(t\right)exp(S^{★}+W^{★}-I^{★})-(\alpha +h+g)\left(t\right),\hfill \end{array}$$

(62b)

$$\begin{array}{cc}\hfill & \frac{d{U}^{★}}{dt}=a\left(t\right)(1+exp(W^{★}-S^{★}))\left(1-\frac{exp\left(U^{★}\right)+exp\left(W^{★}\right)}{K(exp\left(S^{★}\right)+exp\left(I^{★}\right))}\right)\hfill \\ \hfill & -\gamma \left(t\right)exp\left(I^{★}\right)-(\mu +c)\left(t\right),\hfill \end{array}$$

(62c)

$$\begin{array}{cc}\hfill & \frac{d{W}^{★}}{dt}=\gamma \left(t\right)exp(I^{★}+U^{★}-W^{★})-(\mu +c)\left(t\right),\hfill \end{array}$$

(62d)

$$\begin{array}{cc}\hfill & (S^{★},I^{★},U^{★},W^{★})\left(0\right)=(ln\left(S_0\right),ln\left(I_0\right),ln\left(U_0\right),ln\left(W_0\right)).\hfill \end{array}$$

(62e)

In particular, it follows that the following two statements are valid:

(a) 

If $(S^{★},I^{★},U^{★},W^{★}),$ is ω-periodic, then $(S,I,U,W)$ is ω-periodic.

(b) 

If (62) has a solution, then (13)–(17) has a positive solution.

Proof. 

We prove that systems (13)–(17) and (62) are equivalent, in the sense that $(S,I,U,W)$ is a solution of system (13)–(17) if and only if $(S^{★},I^{★},U^{★},W^{★})$ is a solution of system (62). Firstly, let us consider that $(S,I,U,W)$ is a solution of system (13)–(17) and prove that $(S^{★},I^{★},U^{★},W^{★})$ is a solution of system (62). We observe that multiplying (13) by $exp(-S^{★}\left(t\right))$ and using the change of variable (61), we deduce that

$$\begin{array}{cc}\hfill \frac{d{S}^{★}}{dt}& =\frac{dS}{dt}exp(-S^{★}\left(t\right))\hfill \\ \hfill & =[b\left(t\right)r_S(exp\left(S^{★}\left(t\right)\right),exp\left(I^{★}\left(t\right)\right))exp\left(S^{★}\left(t\right)\right)\left(1-\frac{exp\left(S^{★}\left(t\right)\right)+exp\left(I^{★}\left(t\right)\right)}{\theta }\right)\hfill \\ \hfill & -\beta \left(t\right)exp\left(S^{★}\left(t\right)\right)exp\left(W^{★}\left(t\right)\right)-h\left(t\right)exp\left(S^{★}\left(t\right)\right)+g\left(t\right)exp\left(I^{★}\left(t\right)\right)]exp(-S^{★}\left(t\right))\hfill \\ \hfill & =b\left(t\right)r_S(exp\left(S^{★}\left(t\right)\right),exp\left(I^{★}\left(t\right)\right))\left(1-\frac{exp\left(S^{★}\left(t\right)\right)+exp\left(I^{★}\left(t\right)\right)}{\theta }\right)\hfill \\ \hfill & -\beta \left(t\right)exp\left(W^{★}\left(t\right)\right)-h\left(t\right)+g\left(t\right)exp\left((I^{★}-S^{★})\left(t\right)\right).\hfill \end{array}$$

Then, $S^{★}$ satisfies (62a). Analogously, by using the change of variable defined in (61) and multiplying (14), (15), and (16) by $exp(-I^{★}\left(t\right)),$ $exp(-U^{★}\left(t\right))$, and $exp(-W^{★}\left(t\right)),$ respectively, we deduce that $I^{★},U^{★},$ and $W^{★}$ satisfy (62b), (62d), and (62d). Consequently, $(S^{★},I^{★},U^{★},W^{★})$ is a solution of system (62). Conversely, if we assume that $(S^{★},I^{★},U^{★},W^{★})$ is a solution of (62), we deduce that $(S,I,U,W)$ is a solution of (13)–(17) by (61) and differentiation. For instance, for the first equation, we have that

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =exp\left(S^{★}\left(t\right)\right)\frac{d{S}^{★}}{dt}\hfill \\ \hfill & =exp\left(S^{★}\left(t\right)\right)[b_S(exp\left(S^{★}\right),exp\left(I^{★}\right),t)\left(1-\frac{exp\left(S^{★}\right)+exp\left(I^{★}\right)}{\theta }\right)\hfill \\ \hfill & -\beta \left(t\right)exp\left(W^{★}\right)-h\left(t\right)+g\left(t\right)exp(I^{★}-S^{★})]\hfill \\ \hfill & =b_S(S,I,t)S\left(1-\frac{S+I}{\theta }\right)-\beta \left(t\right)SW-h\left(t\right)S+g\left(t\right)I.\hfill \end{array}$$

To prove (a) and (b), we proceed as follows. If we consider that the functions $S^{★}$, $I^{★}$, $U^{★}$, and $W^{★}$ are $\omega$-periodic, then by (61) we can obtain that the functions $S,$ $I,$ $U,$ and W are $\omega$-periodic, since (for instance, in the case of S) we have that $S(t+\omega )=exp\left(S^{★}(t+\omega )\right)=exp\left(S^{★}\left(t\right)\right)=S\left(t\right)$. Hence, we can prove (a). The item (b) is a straightforward consequence of the change of variable (61). □

4.2. Existence of a Periodic Solution for (62)

Lemma 6. 

Assuming that the coefficients of system (62) satisfy the hypotheses of Theorem 1, then system (62) has at least one ω-periodic solution.

Proof. 

In order to analyze system (62) by applying the topological degree arguments, we reformulate this system (62) as an operator equation. We consider the operators $L:\mathrm{Dom}L\subset X\to Y$ and $N:X\to Y$ as defined by the relations in (21) and (22), with $X=Y$ being the space defined in (20), and $\Omega$ as defined in (23), where $\sigma$ is given in (60). With this setup, system (62) can be rewritten as the operator equation:

$$\begin{array}{c}\hfill L{\left((S^{★},I^{★},U^{★},W^{★})\right)}^T=N{\left((S^{★},I^{★},U^{★},W^{★})\right)}^T,(S^{★},I^{★},U^{★},W^{★})\in \mathrm{Dom}\left(L\right)\subset X.\end{array}$$

(63)

Then, by Lemma 4 and Theorem 2 we deduce that system (62) has at least one solution in $\mathrm{Dom}L\cap \overline{\Omega }\subset X$. □

4.3. Proof of Theorem 1

From Lemma 5, we observe that system (13)–(17) can be rewritten as system (62). As a consequence of Lemma 6, we deduce that system (62) has at least one $\omega$-periodic solution. Therefore, by items (a) and (b) of Lemma 5, we conclude that system (13)–(17) has at least one solution that is both $\omega$-periodic and positive.

5. Numerical Examples

In this section, a series of numerical examples are presented. The model study is motivated by the cassava mosaic virus, which is carried by the whitefly Bemisia tabaci; the simulations have the goal of visualizing general model behavior that goes beyond a specific virus. Therefore, the approach can be applied to other types of crops where a cutting (like a stem, branch, or shoot) is used for grafting. Regarding the units for the parameter values, they vary based on the specific parameters that are listed in Table 1, and where the units are given. For example, the maximum replanting rate is measured in units per day, the maximum plant population pressure is measured per square meter, the contact rate between a healthy plant and an infected vector, the harvesting rate, and the recovery rate from infection are all measured in units per day. The parameter values can be adapted to an a-dimensional variable setting, where the variables are normalized by the constraints

$$\begin{array}{c}\hfill S+I\le 1,U+W\le 1.\end{array}$$

(64)

In this a-dimensional variable setting, the parameters can be manually adjusted to illustrate the qualitative model behavior. The general purpose of this contribution is to provide applicability for the non-dimensional case, that can be then adapted to any specific quantitative setting. Regarding specific parameters, there are plants with similar live cycle as the cassava plant, e.g., sweet potatoes.

5.1. Parameter Adjustment

Replication cycles of the cassava virus are studied in [38,39]. Viruses themselves do not have a “life cycle”, as they require a host to replicate. In the whitefly life cycle, the adulthood spans from their second to fourth week, leading to multiple generations per year. The baseline life cycle though is that of the cassava plant: the virus can persist in an infected plant for its entire life span, which is typically 1–3 years, depending on the variety and cultivation practices. The literature tells us that, after planting, the cassava harvest cycle is between the fourth and ninth month [40]. So, as a mean, a reference period of 200 days is chosen.

To give some orientation of the numerical parameters, we compare the parameters of our model in Equations (13)–(17) with the parameters in some studies in the literature, namely, those of Bahity-Maranya (2022) [41] (model 2),

$$\begin{array}{cc}\hfill \frac{d{S}_c}{dt}& =r_c\left(1-\frac{N_c}{{\kappa }_c}\right)S_c+\varphi I_c-a{S}_c{I}_w-\psi S_c\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill \frac{d{E}_c}{dt}& =r_c\left(1-\frac{N_c}{{\kappa }_c}\right)\rho E_c+a{S}_c{I}_w-(\psi +ϵ)E_c\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill \frac{d{I}_c}{dt}& =ϵE_c-(\psi +\varphi )I_c\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill \frac{d{S}_w}{dt}& =\left(1-\frac{N_w}{{\kappa }_w}\right)r_w{N}_w-\left(b(E_c+I_c)+c{I}_h\right)S_w-\omega S_w\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill \frac{d{I}_w}{dt}& =\left(b(E_c+I_c)+c{I}_h\right)S_w-\omega I_w\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill \frac{d{S}_h}{dt}& =r_h\left(1-\frac{N_h}{{\kappa }_h}\right)S_h-d{S}_h{I}_w-\eta S_h\hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill \frac{d{I}_h}{dt}& =d{S}_h{I}_w-\eta I_h\hfill \end{array}$$

(71)

and Magoyo et al. (2019) [42] (model 3),

$$\begin{array}{cc}\hfill \frac{d{S}_r}{dt}& =r_1{S}_r\left(1-\frac{S_r}{k_1}\right)-{\beta }_1{S}_r{I}_V-{\rho }_1{S}_r\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill \frac{d{S}_C}{dt}& =r_2{S}_C\left(1-\frac{S_C}{k_2}\right)-{\beta }_2{S}_C{I}_V-{\rho }_2{S}_C\hfill \end{array}$$

(73)

$$\begin{array}{cc}\hfill \frac{d{I}_C}{dt}& ={\beta }_2{S}_C{I}_V+{\beta }_1{S}_r{I}_V-{\rho }_3{I}_C-a{I}_C\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill \frac{d{S}_V}{dt}& =b(S_V+I_V)\left(1-\frac{S_V+I_V}{k_3}\right)-{\beta }_3{S}_V{I}_C-\gamma S_V\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill \frac{d{I}_V}{dt}& ={\beta }_3{S}_V{I}_C-\gamma I_V\hfill \end{array}$$

(76)

The parameters of these models are compared in the concordance in

Table 2. Concordance table comparing parameters of three models.

	Model 1			Model 2 [41]			Model 3 [42]
	(System (13)–(17))			(System (66)–(71))			(System (72)–(76) )
$b\left(t\right)$	$0.05{\mathrm{day}}^{-1}$		$r_c$	$0.05{\mathrm{day}}^{-1}$		$r_1$	$0.025{\mathrm{day}}^{-1}$
				$0.2{\mathrm{day}}^{-1}$		$r_2$	$0.2{\mathrm{day}}^{-1}$
$\theta$	$1.0{\mathrm{plant}/\mathrm{m}}^2$		${\kappa }_c$	$0.7{\mathrm{m}}^{-2}$		$k_1$	3000
						$k_2$	2000
$\beta \left(t\right)$	$0.001{\mathrm{vector}}^{-1}{\mathrm{day}}^{-1}$		a	$0.008{\mathrm{plant}}^{-1}{\mathrm{day}}^{-1}$		${\beta }_1$	$0.0012{\mathrm{vector}}^{-1}{\mathrm{day}}^{-1}$
						${\beta }_2$	$0.003{\mathrm{vector}}^{-1}{\mathrm{day}}^{-1}$
$h\left(t\right)$	$0.03{\mathrm{day}}^{-1}$		$\psi$	$0.003{\mathrm{day}}^{-1}$		${\rho }_1$	$0.005{\mathrm{day}}^{-1}$
						${\rho }_2$	$0.003{\mathrm{day}}^{-1}$
$g\left(t\right)$	$0{\mathrm{day}}^{-1}$
$a\left(t\right)$	$0.1{\mathrm{day}}^{-1}$		$r_h$	$0.02{\mathrm{plant}}^{-1}{\mathrm{day}}^{-1}$		b	$0.5{\mathrm{vector}}^{-1}{\mathrm{day}}^{-1}$
$\gamma \left(t\right)$	$0.01{\mathrm{plant}}^{-1}{\mathrm{day}}^{-1}$		d	$0.008{\mathrm{day}}^{-1}$		${\beta }_3$	$0.002{\mathrm{plant}}^{-1}{\mathrm{day}}^{-1}$
$\mu$	$0.02{\mathrm{day}}^{-1}$		$\eta$	$0.001{\mathrm{day}}^{-1}$		$\gamma$	$0.0782{\mathrm{day}}^{-1}$
$\alpha$	-

, where corresponding parameters are placed in the same row.

5.2. Specific Examples

The examples are divided based on their underlying scenarios: examples 1 and 2 demonstrate the fundamental behavior in non-periodic situations, while examples 3 and 4 focus on periodic situations. For all the examples presented here, we have consistently used the following initial condition:

$$(S_0,I_0,U_0,W_0)=(0.1,0.1,0.1,0.1).$$

5.2.1. Example 1

In the first example, we start with standard values that a oriented at the literature, see Table 2, and compare variations in abundance replanting with frequency replanting. Specifically, different variants of abundance replanting and frequency replanting are simulated, namely, abundance replanting with $\widehat{ϵ}\in \{0.5,0.8\}$ (-), and frequency replanting with $\epsilon \in \{0.7,0.8\}$ (-).

Figure 2 illustrates that all cases exhibit similar qualitative behavior: The sum of S and I approaches a carrying capacity. The value of S initially increases, reaching a maximum in cases 1a and 1c, whereas in case 1b, the decrease is only slight, and in 1d, S appears to increase monotonically within the visualized time window. The variable I behaves in a manner complementary to S: it increases monotonically in 1a and 1c. In 1b, it rises to a local maximum, then slightly decreases to a local minimum, before increasing slightly again.

The situation for U is clearer; in all examples, it increases to a maximum before subsequently decreasing. Complementarily, W exhibits S-shaped curves that resemble the behavior of a logistic equation.

In subsequent examples, specific parameters are varied, while the other variables are kept at reference levels.

5.2.2. Example 2

In example 2, we use the parameter set from example 1a as a baseline, which corresponds to the case of abundance replanting with $\widehat{ϵ}=1$ (-).

From this baseline, each example involves altering a single parameter as follows: $\beta =0.01$ (vector${}^{-1}$ day${}^{-1}$), $g=0.01$ (day${}^{-1}$), $h=0.01$ (day${}^{-1}$), and $\mu =0.01$ (day${}^{-1}$), as depicted in Figure 3.

The increase in $\beta$ to 0.01 (vector${}^{-1}$ day${}^{-1}$) results in I dominating S early on, yet both variables grow in parallel, as illustrated in Figure 3a.

The higher recovery rate $g=0.01$ (day${}^{-1}$), which models recovery from infection, not only tempers the rise in I but also allows for strong early growth in S. This enables U to reach and maintain a relatively high asymptotic value after surpassing its peak, as shown in Figure 3b.

Implementing continuous harvesting by setting $h=0.01$ (day${}^{-1}$) moderates the increases in both S and I, as shown in Figure 3c. Notably, higher values such as $h=1$ (day${}^{-1}$) would visibly constrain and prolong their development.

Finally, increasing vector species’ mortality to $\mu =0.01$ (day${}^{-1}$) allows U to reach higher values, giving more space to S and moderating the level of I, as seen in Figure 3d.

5.2.3. Example 3

The parameters in example 3 keep the baseline values of example 1a, with the only modification being a periodic harvesting parameter defined as

$$\begin{array}{c}\hfill h\left(t\right)=\left\{\begin{array}{cc}0{\mathrm{day}}^{-1}\hfill & \mathrm{if}tmod250\le 240\mathrm{days},\hfill \\ 0.1{\mathrm{day}}^{-1}\hfill & \mathrm{if}tmod250>240\mathrm{days}.\hfill \end{array}\right.\end{array}$$

(77)

We refer to Figure 4 for illustration: The harvesting leads to a simultaneous reduction in S and I, which is accompanied by a reduction in U and W. The decrease in W towards a minimum allows S to obtain a maximum, that is accompanied by a maximum in U. After reaching their maxima, S and U decrease towards an asymptotic state, whereas I and W continue to increase towards an asymptotic state.

Increasing the harvesting parameter changes the phenomenon quantitatively, but not qualitatively. Remarkable, higher harvesting parameters $h\left(t\right)$ might lead to negative U values. Similarly, changing the mortality of the vector species might also result in negative values of U. As negative values of U are not physically plausible, a model modification is necessary to exclude negative values by design. Since the negative values can be traced back to the term in the large parenthesis, there are two possible fixes:

$$\begin{array}{c}\hfill U\left(1-\frac{U+W}{K(S+I)}\right),max(0,\left(1-\frac{U+W}{K(S+I)}\right)).\end{array}$$

(78)

Multiplying the potentially negative factor by U ensures that the values of U remain within the non-negative invariance region, with $U=0$ as an asymptotic boundary that cannot be crossed. Taking the maximum of the potential negative value and 0 enforces a non-negative value.

5.2.4. Example 4

In example 4, we combine the periodic harvesting with the settings of 1d. There, for $r_S(t,S,I)$ and $r_I(t,S,I)$ frequency replanting is considered with $\epsilon =0.5$ (-). The effect is quantitative, leading to lower values of W and I, which, in turn, correspond to higher values of U and S; see Figure 5.

6. Conclusions and Future Work

In this paper, we have presented a unified mathematical model for propagation dynamics of a pathogen living in a vector population that affects a plant population. The model considers the density distribution of both plants and vectors. Specifically, the plant population is categorized into two groups, healthy plants and diseased plants, with their time-dependent density populations denoted by S and I, respectively. Similarly, the vector population is categorized into healthy and diseased vectors, with densities denoted by U and W, respectively. These categorizations are part of the approach to model the dynamics of plant virus disease epidemics that affect plant populations and are transmitted by vector populations.

We have successfully demonstrated the existence of positive periodic solutions, assuming that the coefficients are positive functions and the initial condition is positive. Significantly, our main result constitutes an advancement over the findings previously reported in [35]. The methodology employed to prove this result is based on the coincidence degree theory.

Additionally, we have illustrated the model’s implications through five numerical examples. These examples effectively demonstrate the periodic behavior of the system’s solutions under the conditions stipulated by our main theorem. This practical application highlights the relevance and applicability of our theoretical findings in understanding the dynamics of pathogen propagation in vector–plant interactions.

The real world applications of model predictions supports a predictive agriculture: With better models, farmers can predict when their crops are most at risk and take preventive measures, which leads to more stable food supplies and economic stability for farming communities. Complementarily, such predictions allow for resource allocation, e.g., of pesticides and labor, and biodiversity conservation by predicting and preventing outbreaks. Agricultural disease management and ecosystem health, an understanding of the periodic solutions of plant virus systems, helps in predicting disease outbreaks, optimizing control strategies, and enhancing crop protection.

In our future work, we plan to conduct research on the epidemiology of plant–insect systems in at least three ways. First, we plan to study the comparison of the model with experimental data via the analysis of the parameter identification problem by using steepest descent methods with the adjoint and sensitivity methodologies [43,44]. The second idea is to improve the model by including the effects of diffusion in the vector displacement and using the ideas of the SIS diffusion–reaction model [45]. Moreover, we expect to develop the model in the context of fractal-fractional differential theory [8].