By including this overview in the biological foundations, readers gain an integrated view of disease dynamics and intervention strategies, while the technical and quantitative aspects of the model are addressed later .

1. The causal agent of Citrus Greening: Candidatus Liberibacter asiaticus

The pathogen of Greening is a gram‐negative and phloem‐limited bacterium, and there are three species, Candidatus Liberibacter asiaticus (CLas), Candidatus Liberibacter africanus (CLaf) and Candidatus Liberibacter americanus (CLam), but the major causal agent is CLas^[1,4]. It is the most feared disease among orange producers because there is no cure or treatment and it spreads very quickly ^[5].

Currently, CLas cannot be cultured in vitro due to its unclear growth factors, making it a challenge studying Greening ^[4]. Also, the characteristic of phloem-limited pathogensmakes it especially difficult to control ^[4]. The bacterium lives in the ‘veins’ of the plants (phloem vessels), through which it spreads rapidly to all parts of the tree: roots, branches, leaves and fruits. After the bacterium is transmitted by the psyllid, the symptoms of Greening begin to appear on the tree's leaves between 4 and 10 months. In young plants, the tree can be fully compromised in one or two years, while older plants can take three to five years to be fully affected by the symptoms ^[9] . By colonizing this tissue, the bacterium disrupts the flow of sugars and essential nutrients for plant growth and storage, causing severe systemic symptoms ^[5,6]. Infection leads to a gradual loss of citrus plant productivity, causing yellowing and blotchy mottling, deformed fruits, and reduced nutrient transport^[5,6].

2. Phloem structure and function

The phloem is responsible for interconnecting the whole plant level by long distance two way dual function transport ^[7]. The water, photosynthetic products, proteins and other molecules flow in the phloem driven by an osmotically generated pressure gradient ^[7].

Transport occurs through sieve tube elements, connected by sieve pores that allow substances to pass and ensure continuity of flow within the plant ^[8]

The transport phloem is dual function: distributes assimilates both to terminal sinks (roots, shoots) and to lateral sinks(e.g., cambium) that sustain growth and tissue maintenance. Along the pathway, solutes continuously leak and are retrieved ^[7]. Loading and unloading can be active or passive, occurring either apoplastically (via membrane transport into the apoplast) or symplastically (through plasmodesmata connections) ^[7].

Colonization of the phloem by Clas compromises this transport, leading to sugar accumulation in the leaves and deficiency in terminal organs, contributing to the classic symptoms of Greening ^[2,3].

3. Plant immune response

Upon detecting CLas, the plant activates local defenses in the phloem:

1. Callose deposition: polymer deposited at sieve pores that limits bacterial movement ^[3].
2. Localized immune suppression: CLas can interfere with immune activation, allowing its spread.

These responses form a delicate balance: they attempt to contain the bacterium but can harm the plant itself by blocking nutrient transport or causing cellular damage ^[3].

4. Pathogen strategies

In citrus Greening, the CLas achieves systemic movement by manipulating the plant’s defense responses in the phloem. In that way, CLas overcomes the plant’s physical defenses in the phloem by reducing or bypassing callose accumulation, allowing it to continue spreading systemically despite the plant’s attempt to block its movement ^[2,3].

Under normal conditions, healthy plants deposit callose around sieve pores as a physical barrier to restrict pathogen spread, a process triggered by calcium signaling, reactive oxygen species (ROS), and salicylic acid^[2,3].

However, CLas suppresses callose deposition and reduces ROS accumulation in infected sieve elements, thereby preventing occlusion of sieve pores and allowing unhindered bacterial passage through the phloem network. This suppression involves downregulation of callose synthase (CalS) and phloem protein 2 (PP2) genes, as well as the secretion of bacterial effectors such as peroxiredoxin, which mitigate ROS toxicity and dampen immune signaling. ^[2,3]

Thus, the pathogen’s mobility is tightly linked to its ability to disrupt the ROS–callose defense axis, a strategy that shifts the balance from host containment toward bacterial dissemination. In that way, while strong induction of these defenses would normally confine pathogens, their inhibition by CLas ensures both survival and systemic colonization of the phloem ^[2,3].

5. Therapeutic Intervention Module: Rationale and Simulation Strategy

Our model was designed to evaluate and compare therapeutic intervention strategies against Greening. To achieve this, we simulated the action of two compounds with distinct mechanisms of action: our candidate peptide CTX and oxytetracycline ^[11], a benchmark antibiotic used in the field, specially in the United States.

The Candidate Peptide (CTX) a Bactericidal Action

The primary strategy is based on the action of antimicrobial peptides (AMPs), a class of molecules recognized in the literature for their potent activity against a wide range of plant pathogens, including difficult-to-control bacteria ^[12]. The predominant mechanism of action for many AMPs is bactericidal, resulting in the direct death of the pathogen through membrane destabilization or the disruption of vital metabolic processes^[13].

Based on this biological rationale, our model assumes that CTX operates through a similar mechanism. Mathematically, this is implemented in Module 3 as a mortality term added to the population dynamics equation of CLas, actively removing a fraction of the bacterial population at each time step.

The Benchmark Treatment (Oxytetracycline) – Bacteriostatic Action

To validate our model and establish a realistic benchmark, we also simulated the action of oxytetracycline. This simulation was parameterized based on the commercial product ArborBiotic™ by INVAIO, which uses oxytetracycline for Greening control. According to the product documentation, its mode of action is the inhibition of protein synthesis, which characterizes a bacteriostatic action; it prevents bacterial replication but does not cause direct cell death ^[11].

To anchor our simulation in a real world agronomic scenario, we adopted the maximum recommended dosage for mature trees: 150 mg of oxytetracycline per season. In the model, the bacteriostatic action is implemented differently from the bactericidal one: instead of adding a mortality term, we reduce the intrinsic growth rate in the logistic equation of the CLas population.

Pharmacodynamic Implementation and Model Applications

For both treatments, the relationship between the compound's concentration and its effect (whether bactericidal or bacteriostatic) is modeled using the Hill function, a standard and widely validated mathematical approach in pharmacology for describing saturable biological responses ^[14]. Furthermore, the dynamics of each compound's concentration over time are governed by pharmacokinetic (PK) principles to simulate realistic exposure profiles ^[15].

This dual approach, grounded in both the general literature on AMPs and specific data from a commercial product like INVAIO's ^[11], allows the model to examine how different interventions can mitigate infection, restore nutrient transport and prevent systemic spread. It becomes a robust tool for comparing the relative efficacy of novel strategies, such as CTX, against an established benchmark treatment. Future applications could extend to other phloem-inhabiting diseases, using specific data to predict disease progression and support the development of targeted treatments.

1. Overview

To investigate the dynamics of Candidatus Liberibacter asiaticus (CLas) in the phloem and evaluate our intervention strategy with CTX, we developed a hybrid cellular automata model ^[1]. This approach represents the phloem as a two-dimensional grid (lattice) of sites with periodic boundary conditions, where each site corresponds to a sieve tube element and its local sieve pore connections.

This abstraction focuses on the transport dynamics relevant to CLas infection and callose deposition, rather than the complete physiology of companion cells or other phloem components. For clarity, throughout the text we use the term “site” to denote these representative phloem elements.

Why a 2D Grid? A Strategic Abstraction of the Phloem

While the phloem is a complex three dimensional network, we made a deliberate design choice to represent it as a 2D grid ^[2]. This is not a literal depiction but a functional abstraction designed to isolate the core mechanisms of disease progression. One can imagine our grid as a "functional plan view"of the vascular tissue, where the most critical interactions—the local, cell-to-cell spread of the pathogen are explicitly modeled.

This strategic simplification provides three advantages:

• Focus on the Essential Mechanism: The spread of CLas is fundamentally a process of local contagion, moving from one sieve tube element to its immediate neighbors through sieve pores. Our 2D grid, with its defined neighborhood structure, captures this core dynamic effectively, allowing us to study the interplay between bacterial propagation and the plant's localized immune response without the confounding complexities of 3D vascular architecture.
• Enabling Robust Analysis: A computationally efficient 2D model allows for extensive simulation runs. This was critical for performing a thorough sensitivity analysis, where we tested how our conclusions hold up under a wide range of parameter values. This ensures that our findings on the efficacy of CTX are robust and not merely an artifact of a single set of assumptions.
• Clarity and Interpretability: The 2D format facilitates the visualization and interpretation of emerging spatial patterns, such as infection waves, callose blockages and zones of recovery. This makes the model a tool not only for quantitative analysis but also for communicating our findings to a broader audience.

While this approach is powerful for studying local dynamics, we acknowledge that it does not capture long distance transport phenomena (e.g., source-sink dynamics) within the plant. This represents a clear avenue for future model enhancements.

Periodic Boundary Conditions: Simulating a Continuous Tissue

We implemented periodic boundary conditions in our grid. This means that the left and right edges, as well as the and bottom edges, are connected. When a process, like bacterial spread, exits one side of the grid, it re-enters on the opposite side. This design choice is critical for two reasons:

1. Eliminating Edge Effects:In a non-periodic grid, sites at the edges have fewer neighbors, creating artificial boundaries that don't exist in a continuous biological tissue. These boundaries can unnaturally halt or alter the propagation patterns of the infection.
2. Approximating an Infinite System: Our simulation represents a small patch of a much larger phloem network. By connecting the edges, we simulate the behavior of this patch as if it were part of an effectively infinite, continuous system, making our results more generalizable and less dependent on the specific size of our grid.

To capture the system's complexity, our model is hybrid:

1. Discrete states to describe the general status of the site (susceptible, infected, or responding with callose);
2. Continuous variables to quantify bacterial concentration and callose deposition.

The model's dimensions allow for the simulation of several key processes:

1. CLas propagation in the phloem;
2. Plant immune responses like callose
3. Therapeutic interventions, such as CTX application.

Ultimately, this combination enables the model to represent gradual transitions and heterogeneous local responses, providing a tool to explore complex plant-pathogen interactions.

2. System Construction

The system is built on a network of cellular automata, where each site in the grid represents a phloem sieve tube element. Local neighborhoods reflect the natural connections of sieve pores, through which bacteria, nutrients, and signaling molecules circulate.

A key feature of our model is its dual-scale nature, reflecting two distinct biological ranges of interaction. Bacterial infection operates on a local, cell-to-cell scale, governed by the immediate Von Neumann neighborhood (left, right, top, and bottom). In contrast, the plant's immune response (callose deposition) is activated on a broader scale, responding to a signal that propagates up to a fixed radius around infected sites. This design realistically captures how defense mechanisms can be triggered in nearby cells even before the pathogen arrives.

2.1. Discrete States

Each site can be in one of three discrete states:

• S - Susceptible;
• I - Infected;
• C - Callose Response;

This classification is not a stored variable but an interpretation of the underlying continuous dynamics, similar to how epidemiological models like SIR categorize individuals ^[3].

2.2. Continuous Variables

In addition to the discrete states, each site carries continuous variables that evolve over time:

• $I_{ij}$ - CLas infection concentration;
• $C_{ij}$ - callose amount;

These variables allow gradual representation of infection intensity and defense responses, instead of only abrupt state changes.

They allow classification of sites into the discrete categories S, I and C based on dynamics. In practice, this means that callose deposition is graded, but for interpretation we describe sites as “C” when the defense response is strong enough to significantly reduce permeability and hinder infection. Otherwise, the infection spreads and the classification of the site is “I” or “ S”if it is protected by callose. Similar to epidemiological models as SIR ^[3].

2.3. Initial conditions

All sites start in the S state with continuous variable values set to zero ($I_{ij} = 0$ and $C_{ij} = 0$). One or more sites receive an initial concentration of $I_{ij}(0) > 0$, representing the CLas inoculation point.

2.4 Neighborhood and Local Interactions

• Infection dynamics are governed by the Von Neumann neighborhood (up, down, left, right)^[1,2]. An infected site can transmit bacteria to any of its four immediate neighbors.
• Callose activation responds not only to direct infection but also to a diffusive signal extending up to a fixed radius $S_r$(e.g. 2 sites around an infection point). This reflects that phloem defense signaling isn’t restricted to the directly infected site but propagates through short-distance signals, activating callose in neighboring regions before the bacteria arrive ^[4].

2.4. State Transitions

Although implemented with continuous dynamics, the categories evolve conceptually as:

• Infected (I): if $I_{ij} > 0$ and local defenses are insufficient.
• Callose Response (C): if $C_{ij}$ approaches saturation ($C_\text{limit}$).
• Susceptible (S): otherwise, no infection and no defense activation.

All sites are updated simultaneously at each time step ($t \rightarrow t+1$), ensuring temporal consistency and avoiding propagation bias, as infection and defense activation occur in parallel.

3. Auxiliary Functions

In our model, we incorporate auxiliary functions that represent how the CLas spreads through the phloem and how the plant activates its defenses. These functions describe concentration dependent responses, supported by experimental evidence. For example, in sweet orange trees infected by CLas defenses aren’t activated immediately after infection but only when the bacterial concentration exceeds critical levels ^[5]. This delay results in late activation followed by physiological repression, a pattern also reported in other studies regarding callose in the phloem ^[6,7].

This behavior can be described by nonlinear mathematical functions, such as the Hill function (for defense activation), and an exponential function (for permeability reduction by callose accumulation that regulates the probability of infection). These functions capture gradual and saturable transitions rather than abrupt changes. If you want to see more about these functions and some examples, you can access our supplementary material.

3.1. Modeling Threshold-Based Responses: The Hill Function

Many biological processes, such as immune activation or drug efficacy, do not respond linearly to stimuli. Instead, they often exhibit a switch-like behavior, where the response is minimal until a critical concentration (a threshold) is surpassed, after which it activates rapidly. To mathematically capture this sigmoidal behavior, our model uses the Hill function. This versatile function is applied in two contexts ^[8,9]:

Module 2 (Plant Response): To model the activation of callose deposition in response to the local bacterial signal.

Module 3 (Therapeutic Interventions): To model the pharmacodynamic relationship between drug concentration and its effect (either bactericidal or bacteriostatic).

The general form of the Hill function $H(s)$ used in our model is:

\begin{equation} H(s) = \frac{s^n}{s_{50}^n + s^n} \end{equation}

Where:

• $s$: Represents the concentration of the input signal (e.g., the local bacterial concentration $I_{\text{local}}$ for callose, or the drug concentration $D(t)$ for therapies).
• $s_{50}$: The activation constant (also known as EC$_{50}$), representing the signal concentration at which the response reaches 50% of its maximum. This is the activation threshold.
• $n$: The Hill coefficient, which defines the steepness and sensitivity of the response. Higher values of $n$ result in a more switch-like behavior.

This function returns a value between 0 (no response) and 1 (maximum response), effectively representing how a biological system can shift from an "off" to an "on" state once a stimulus surpasses a critical level.

3.2 Sieve Pore Permeability: Exponential Function

The primary effect of callose in the phloem is to restrict the permeability of sieve pores, limiting the transport of both nutrients and pathogens ^[10,11]. We represent this process mathematically as an exponential decay function:

\begin{equation} P_{ij}(t) = e^{-\lambda_C \cdot C_{ij}(t)} \end{equation}

Where:

• $C_{ij}(t)$: amount of callose in site $(i,j)$;
• $\lambda_C$: sensitivity of permeability to callose accumulation.

This function ensures that even at high callose levels, permeability is never fully zero but decreases sharply as deposition increases. Biologically, this reflects the fact that callose can strongly hinder CLas propagation without necessarily achieving complete blockage of the spread of the disease.

In our model, we use this permeability function directly regulates the probability of infection:

\begin{equation} P_{\text{infection}}^{ij} \propto P_{ij}(t) \end{equation}

Thus, higher local callose concentrations reduce the success of bacterial colonization. This mechanistic link grounds the model in experimental evidence showing that callose accumulation in sieve pores is a major factor controlling pathogen spread in the phloem ^[11].

The functional forms of these relationships are illustrated conceptually below. Panel (A) shows the exponential decay in permeability as callose increases, corresponding to Equation 2. Panel (B) shows the direct linear proportionality between permeability and the resulting infection probability, as described in Equation 3.

4. Modules

4.1 Module 1: CLas Bacterial Propagation

In our model, the propagation of CLas in the phloem is a stochastic discrete spreading to neighboring sites combined with continuous growth of the bacterial population at each site.

We determine the bacterial concentration $I_{ij}(t)$ in each site as:

\begin{equation} I_{ij}(t+1) = I_{ij}(t) + I_{ij}(t) r \left( 1 - \frac{I_{ij}(t)}{I_{\max}} \right) - d C_{ij}(t) I_{ij}(t) - \delta_I I_{ij}(t) \end{equation}

Where:

• $I_{ij}(t)$: bacterial load (infection) at site $(i,j)$ at time $t$;
• $r$: intrinsic bacterial growth rate;
• $I_\text{max}$: maximum bacterial capacity per site;
• $C_{ij}(t)$: callose concentration at site $(i,j)$;
• $d$: suppression coefficient representing the inhibitory effect of callose;
• $\delta_I$ rate of natural degradation;

Each term of this equation has a biological interpretation. We based our fundamentation consider some topics:

1. Continuous logistic growth: $r I_{ij} \left( 1 - \frac{I_{ij}}{I_\text{max}} \right)$

This term represents the bacterial replication at each site, which corresponds to a sieve tube element. Logistic growth limits the bacterial population within each site.If there were no resource limitations or callose, the population would grow exponentially. Therefore, $I_{max}$ is the carrying capacity (maximum saturation of the bacterial population in a site).

2. Callose suppression: $-dC_{ij} I_{ij}$

Callose mortality or suppression term. This tells us that host defense reduces the bacterial population proportionally to the amount of callose present.

3. Natural degradation: $-\delta_I I_{ij}$

This term represents the natural degradation or death rate of bacteria. It indicates that a constant fraction of the bacterial population ($I_{ij}$​) is removed at each time step. Biologically, this can be interpreted as the end of the bacteria's lifespan or their passive clearance by the plant's phloem sap flow.

4.1.1 Discrete Spread

Bacteria can spread from an infected site to neighboring sites with low callose concentration.

The probability of infection of a neighboring site $(i,j)$ is given by:

\begin{equation} P_{\text{infection}}^{ij}(t) = \beta P_{ij}(t) = \beta e^{-\lambda_C C_{ij}(t)} \end{equation}

Where $\beta$ is the base propagation rate and $C_{ij}(t)$ is the local callose concentration.

Each infected site attempts to spread to its neighboring sites independently. For each neighbor that is not yet infected, a stochastic trial is performed: a uniform random number $U \sim \mathcal{U}(0,1)$ is drawn, and the neighbor becomes infected only if $U < P_{\text{infection}}^{ij}(t) = \beta \, e^{-λ_C C_{ij}}$. If the trial succeeds, a small infection seed $I_{\text{seed}}$ is placed in the target site. This process is repeated for each infected site and each of its neighbors independently.

As callose increases, $P_{\text{infection}}$ becomes very small, making it unlikely that the random draw succeeds; with low callose, the probability is higher, so the disease spreads more easily.

When the stochastic trial succeeds, a small infection seed $I_{\text{seed}}$ is placed in the newly infected site.

In summary:

• The previous state $I(t)$ is used to determine potential new infections;
• Successful trials establish $I_{\text{seed}}$ in the target site (without summing multiple seeds);
• The continuous dynamics (logistic growth with callose suppression) are then applied to obtain $I(t+1)$.

4.1.2 Propagation Scenarios

Our model allows us to simulate emergent patterns, such as infection waves, zones protected by callose and areas vulnerable to systemic propagation.

To illustrate these mechanisms, we generated spatial patterns of callose deposition with irregular and scattered formations. These formations mimic the way callose sometimes appears in discrete clusters or corner-like blocks that isolate regions of the phloem.

A key factor in these scenarios is the pore permeability ($P_{ij}$): even with the same level of callose, differences in $P_{ij}$ strongly affect how easily bacteria can spread between sites.

• No Callose: fully open network, representing unrestricted bacterial propagation. (If $C_{ij}(t) = 0$ and $P_{ij}(t) = 1$ (maximum permeability), pores remain fully open and bacteria propagate freely between sites.)
• Partial Callose Response: scattered formations that reduce permeability but still leave “leakage” paths, allowing some bacterial spread. (If $0 < C_{ij}(t) < 1$ and $0 < P_{ij}(t) < 1$, callose deposition partially narrows pores, reducing but not fully blocking bacterial entry and expansion.)
• Full Callose Blockage: large scattered patches that can effectively isolate entire regions of sites, acting as strong barriers. (If $C_{ij}(t) \approx 1$ (high callose) or $P_{ij}(t) \ll 1$, pores are nearly sealed, bacterial entry is strongly hindered, and propagation is effectively restricted.)

The visualizations below demonstrate how these configurations influence CLas dynamics in the phloem. While uniform callose layers are rarely observed biologically, scattered and irregular clusters may create local bottlenecks that slow down or redirect pathogen movement. In all cases, pore permeability is a critical determinant of whether local infections can escalate into systemic spread.

We consider three cases: (i) no callose formation, (ii) partial callose deposition with varying pore permeability, and (iii) complete callose blockage with fully impermeable pores.

These scenarios highlight the importance of permeability in determining how far the disease can spread, whether it can pass through different regions, and how callose affects the increase of average bacterial concentration across the network.

Although callose reduces and often hinders pathogen spread, the disease can sometimes still move into multiple branches or, in other cases, remain localized. The outcome depends strongly on the parameter values discussed previously.

4.2 Module 2: Plant Response (Callose Deposition)

Callose is a polysaccharide that accumulates in sieve pores of the phloem, functioning as a physical barrier that restricts the spread of CLas between sites. Unlike the static barriers considered in Module 1, callose deposition is described here as a dynamic immune response. It’s induced not only at infected sitesbut also in neighboring sites that perceive infection-related signals. This transient and spatially regulated behavior requires callose to be modeled as a time dependent process that balances local production, natural degradation and saturation effects.

4.2.1 Effective Local Signal

In our model, a site $(i,j)$ doesn’t respond only to its own infection but also integrates signals from a surrounding neighborhood.

The local bacterial signal perceived by the site is defined as the average bacterial concentration over a neighborhood of radius $R$ (Manhattan distance) and we describe as:

\begin{equation} I_{\text{local}}^{ij}(t) = \frac{1}{|N_R(i,j)|} \sum_{(m,n) \in N_R(i,j)} I_{m,n}(t) \end{equation}

Where:

• $N_R(i,j)$ is the set of all sites $(m,n)$ within Manhattan distance $R$ from $(i,j)$;
• $I_{m,n}(t)$ is the bacterial concentration at site $(m,n)$ at time $t$.

This formulation captures the idea that phloem elements integrate both direct bacterial load and extracellular infection cues.

4.2.2 Dynamics of Callose

We formulate that the dynamics of callose deposition at site $(i,j)$ is described by:

\begin{equation} C_{ij}(t+1) = \min \Big[ C_{\max}, \; C_{ij}(t) + \alpha_C \, H(I_{\text{local}}^{ij}(t)) - \delta_C \, C_{ij}(t) \Big] \end{equation}

Where:

• $C_{ij}(t)$: amount of callose at site $(i,j)$ at time $t$;
• $\alpha_C$: maximum rate of callose production;
• $\delta_C$: natural degradation rate of callose;
• $C_{max}$: maximum achievable callose (saturation);
• $H(I_{\text{local}}^{ij}(t))$: Hill function describing activation by bacterial signal.

All terms in this equation can be interpreted biologically, and our approach is based on certain foundational topics

1. Production: $\alpha_C H(I_{\text{local}}^{ij}(t))$

This term increases callose proportionally to the local bacterial signal. The Hill function ensures a saturating response, approaching a maximum rate defined by $\alpha_C$.

2. Degradation: $-\delta_C C_{ij}$

This term represents natural turnover of callose over time, reflecting that the plant removes callose to maintain phloem functionality and avoid permanent blockage.

3. Saturation:$\min(\cdot, C_{max})$

This term enforces a biological upper limit, ensuring that callose does not exceed the maximum physically possible deposition $C_{max}$.

4.2.3 Connection with Bacterial Spread

Deposited callose modulates the permeability of sieve pores, which in turn affects bacterial propagation (Module 1). Higher local callose $C_{ij}$ leads to lower pore permeability $P_{ij}$, decreasing the probability of bacterial spread to this site. Conversely, lower $C_{ij}$ leads to higher $P_{ij}$, allowing easier propagation.

This modulation represents the primary interaction between the plant defense response and bacterial growth. By dynamically adjusting permeability according to local bacterial signals, our model captures how spatially heterogeneous callose depositioncan create barriers, slow down infection waves and redirect pathogen spread. In this section demonstrating how local immune responses influence global infection dynamics.

Emergent Behavior

As mentioned before, callose is not a static feature; here we demonstrate the dynamics of its deposition. Our formulation enables the model to reproduce observed characteristics of callose in Greening infected phloem:

• Scattered callose deposition: localized activation produces heterogeneous callose patches.
• Transient blockages: callose accumulates near infection sites but later degrades, reopening pores.
• Incomplete barriers: deposition rarely reaches full occlusion, but pore permeability is significantly reduced, slowing bacterial spread.

4.3 Module 3: Therapeutic Interventions

To validate our therapeutic strategy, our model integrates the action of medicinal agents, allowing us to simulate and compare the effects of different treatments. We adopted two approaches: the bactericidal action of CTX peptide and the bacteriostatic action of tetracycline, an antibiotic sometimes used for controlling Greening.

The distinction between these two mechanisms is fundamental. A bactericidal substance like CTX directly kills bacteria, while a bacteriostatic one like oxytetracycline inhibits their growth. Both actions are modeled in our simulation to capture their respective pharmacodynamics.

The therapeutic agent at the core of our intervention strategy is CTX, a novel antimicrobial peptide synthesized and tested by our team's experimental division. Preliminary in-planta assays confirmed its bactericidal activity (see the full results on our Engineering page). Insights from our model informed the choice of the initial dose for our first in-planta assay. The preliminary results from this assay were inconclusive, and this experimental data is now being used to refine the model's parameters. This feedback loop—from model to experiment and back to the model—is a central component of our project methodology.

For our oxytetracycline treatment simulation, we modeled a regime that reflects a specific agricultural practice. Based on the guidelines of ArborBiotic™ ^[12], a commercial oxytetracycline product used for controlling Greening (Invaio, n.d.), our model simulates an annual single-application protocol recommended for mature orange trees (trunk diameter > 6 cm), which consists of a maximum 150 mg dose per season ^[12]. Adopting this specific scenario ensures our comparison is directly benchmarked against a precise therapeutic strategy used in the field. (See Model Parametrization and Results)

4.3.1 Pharmacokinetics: Drug Concentration

The effectiveness of a treatment depends on its concentration over time. Our model incorporates a simplified pharmacokinetics (PK) module that describes the absorption and elimination of a drug within the plant system. We assume a uniform systemic concentration $D(t)$ that changes over time, governed by parameters such as dose, half-life ($t_{1/2}$), and time to reach maximum concentration ($T_{max}$). This approach allows us to simulate realistic drug exposure profiles ^[13,14].

4.3.2 Pharmacodynamics: Mechanisms of Action

The effect of the drug on the bacterial population is described by adding a new term to the infection dynamics equation from Module 1. The efficacy of the drug is determined by a Hill function, which relates the systemic drug concentration $D(t)$ to its biological activity ^[9].

Bactericidal Action (CTX)

The CTX peptide acts by increasing the bacterial death rate. This is represented by an additional mortality term proportional to the current infection load and the drug's efficacy. The bacterial dynamics equation is modified as follows ^[16]:

\begin{equation} I_{ij}(t+1) = \underbrace{I_{ij}(t) + r I_{ij}(t) \left( 1 - \frac{I_{ij}(t)}{I_{\max}} \right) - d C_{ij}(t) I_{ij}(t) - \delta_I I_{ij}(t)}_{\text{Original Dynamics}} - \underbrace{E_{\text{kill}} \cdot I_{ij}(t)}_{\text{CTX Effect}} \end{equation}

• $E_{\text{kill}}$: The effective killing rate induced by the drug, defined as:

• $D_{\text{CTX}}(t)$: Systemic concentration of CTX at time $t$.
• $EC_{50}$: Drug concentration producing 50% of the maximum effect.
• $n$: The Hill coefficient.
• $\text{killScale}$: The maximum efficacy factor of the bactericidal drug.

Biological Interpretation:

• Drug-Induced Mortality: -$E_{\text{kill}}$​⋅$I_{ij}(t)$ This term represents the direct killing of bacteria by CTX. The rate of killing is proportional to the existing bacterial population Iij and the calculated efficacy $E_{\text{kill}}$, which depends on the drug concentration.

Bacteriostatic Action (Oxytetracycline)

Oxytetracycline acts primarily by inhibiting bacterial growth. In our model, this is captured by reducing the logistic growth term. Additionally, we include an "active clearing" mechanism, representing how the drug facilitates the host's removal of the now-static bacteria ^[12,16].

\begin{equation} I_{ij}(t+1) = I_{ij}(t) + \underbrace{r I_{ij}(t) \left( 1 - \frac{I_{ij}(t)}{I_{\max}} \right) \cdot (1 - E_{\text{inhib}})}_{\text{Inhibited Growth}} - \underbrace{\left( d C_{ij}(t) I_{ij}(t) + \delta_I I_{ij}(t) \right)}_{\text{Natural Death Terms}} - \underbrace{E_{\text{clear}} \cdot I_{ij}(t)}_{\text{Oxytetracycline Clearing}} \end{equation}

Where:

• $E_{\text{inhib}}$: The effective growth inhibition factor, calculated with a Hill function similar to CTX.
• $E_{\text{clear}}$: The active clearing rate, proportional to the inhibition factor: $E_{\text{clear}} = \alpha_{\text{clear}} \cdot E_{\text{inhib}}$
• $D_{\text{Tetra}}(t)$: Systemic concentration of Oxytetracycline.
• $\alpha_{\text{clear}}$: The coefficient for active bacterial clearing.

Biological Interpretation:

• Growth Inhibition: The logistic growth term is multiplied by ($1-E_{\text{inhib}}$). When the drug is effective ( $E_{\text{inhib}\to 1}$ ), bacterial replication is halted.
• Active Clearing: -$E\text{clear}⋅I_{ij}(t)$ This term models the removal of bacteria facilitated by the host system once their growth is inhibited by the drug.

4.3.3 Simulation Scenarios

To investigate the efficacy of our CTX peptide, we compared three treatment scenarios: Control (no treatment), CTX Treatment, and Oxytetracycline Treatment. These scenarios allow us to test our hypotheses and demonstrate how CTX can stand out as an effective alternative for controlling Citrus Greening by comparing its action to a reference antibiotic. Our result visualizations highlight how the infection load and callose response evolve in each of these scenarios.

For our ‘virtual laboratory’ to behave like a real citrus plant, we needed to teach it the rules of biology. In this section, we explain how we used data from scientific literature to define the behaviors of the CLas bacteria and the plant’s defense responses. Each parameter of our model was adjusted to reflect the findings of experimental research, ensuring that our simulation was as realistic as possible.

Our model’s parameterization was guided by a semi-quantitative approach focused on capturing the dynamics and emergent behaviors of the pathogen-host interaction, rather than replicating precise absolute values, many of which are unknown for the CLas-citrus system. The parameters were grouped into three categories:

1. Structural and computational parameters: which define the simulation’s architecture;
2. Biological dynamics parameters: whose values were estimated to ensure biologically plausible relationships.
3. Therapeutic intervention parameters: which were directly informed by laboratory data and commercial product documentation.
4. Constants parameters: which represent the foundational hypotheses and core assumptions of the model.

1. Structural and computational parameters

1. Spatial Framework: To capture the spatial dynamics of infection and defense in the phloem, the model represents the tissue as a 50$\times$ 50 two-dimensional (2D) grid (2500 sites in total), where each site corresponds to an individual cell. The grid size was chosen to balance computational efficiency with the ability to observe the emergence of complex spatial patterns, avoiding finite-size artifacts that could occur in smaller grids. To eliminate edge effects and simulate a portion of a continuous tissue, we implemented periodic boundary conditions.
2. Temporal Dynamics: The simulation progresses in discrete time steps, where each step ($t$) corresponds to one day. The governing equations of the model, therefore, describe the daily dynamics of the system. The total simulation duration is flexible, allowing for the investigation of different time scales, such as the initial phase of infection progression or the long-term effects of a therapeutic intervention. This daily time step is biologically appropriate, as the progression of Greening and the plant's physiological responses, such as callose deposition, occur over days and weeks rather than hours or minutes.
3. Normalization and Numerical Thresholds: To facilitate interpretation and direct comparison between variables, the concentrations of infection (CLas load) and callose were normalized to the [0, 1] interval. On this scale, a value of 1.0 represents the maximum carrying capacity of a site, whether it is the bacterial load the cell can support ($I_{max}$) or the maximum callose deposition ($C_{limit}$). This approach allows other parameters to be defined in relative terms to these maximums.

Parameter Summary

2. Biological dynamics parameters

These parameters govern the core interactions between the Candidatus Liberibacter asiaticus (CLas) bacteria and the citrus host's defense mechanisms. Their values were established through a semi-quantitative approach, designed not to replicate precise numerical values, but to capture the key emergent behaviors of the host-pathogen system as described in the scientific literature.

2.1. Bacterial Dynamics: A Slow and Stealthy Invasion

The progression of Greening is not an aggressive blitz, but a slow, chronic infection characterized by a low bacterial load and an uneven spatial distribution within the host plant [3, 6]. This is supported by findings that CLas is present in only about 7.76% of sieve tube elements in symptomatic leaves^[3], and that the total bacterial population occupies a minuscule fraction, approximately 0.01%, of the total phloem volume in a young tree ^[1]. The CLas genome itself suggests an "obligate parasite" with limited metabolic capabilities, highly dependent on the host for survival and replication ^[3].

Our model captures this behavior through the following parameters:

• Logistic Growth: We use a logistic growth term ($rI(1 − I/I_{max}$​)) to model bacterial replication. This represents a population with limited resources within a host cell, consistent with the long generation time of CLas and its low titer compared to apoplastic pathogens ^[1].
• Propagation Rate ($\beta$): A low base propagation rate was chosen to simulate the slow, localized spread. This is justified by the long latency period observed experimentally, where bacterial populations remain low and stable for up to 90 days post-inoculation before beginning a rapid growth phase ^[4] , and by the consistently low bacterial titers found in young shoots and leaves ^[6].
• Natural Death Rate ($\delta_I$​): A small but non-zero natural death rate reflects the bacterium's nature as an "obligate parasite" that cannot survive independently of the host's cellular machinery ^[3].

2.2. Host Defense Dynamics: A Manipulated and Spatially Displaced Response

The plant's primary defense against phloem-limited pathogens is the deposition of callose, a polymer that functions as a physical barrier in sieve pores ^[5].While Greening is characterized by massive phloem plugging by callose ^[2], recent findings reveal this response is actively manipulated. The defense is spatially displaced: callose levels are low in cells occupied by CLas but high in surrounding, uninfected cells, suggesting that CLas actively suppresses or removes callose to ensure its own mobility ^[3].

Our model integrates these complex dynamics through key parameters:

• Signaling Radius ($S_r$): To reproduce the observed spatial displacement, the signal for callose production affects a neighborhood around an infected cell, matching experimental observations of a strong defensive response in uninfected cells adjacent to infected ones ^[3].
• Callose Production ($\alpha_C$​) and Activation Threshold (Hill Function): The plant's strong "over-reaction" in surrounding tissues ^[2] is modeled using a high maximum production rate ($alpha_C​$) activated by a Hill function. The use of a threshold-based activation is strongly supported by temporal studies showing that disease symptoms, and thus, the strong defense response, only appear after the CLas population reaches a critical concentration of approximately $10^5$ to $10^7$ cells per gram of tissue ^[4]. This aligns with the principles of PAMP-triggered immunity (PTI), which mounts a sharp response when a pathogen signal surpasses a critical threshold ^[5].
• Callose Degradation ($\delta_C$​): Callose deposition is a highly dynamic process involving both synthesis and degradation ^[5]. This parameter models both the natural turnover of callose and, crucially, its active suppression or removal by the bacteria ^[3] . The inclusion of a degradation term is further supported by evidence that the callose response can exhibit a degradation phase after an initial peak, before a massive late-stage accumulation ^[5].
• Callose Efficacy ($d$): To quantify the impact of this physical barrier, this parameter models the capacity of deposited callose to directly suppress bacterial replication by physically blocking sieve pores, as described in the literature ^[5].

Parameter Summary

3. Therapeutic intervention parameters

These parameters define the effects of external treatments designed to combat CLas infection. Our model simulates two distinct therapeutic agents: a custom-designed antimicrobial peptide (CTX) with bactericidal properties [13,14], and a conventional antibiotic, oxytetracycline, known for its bacteriostatic action ^[10]. The parametrization was based on a multiparameter approach to antibiotic treatment modeling ^[11], using established pharmacokinetic and pharmacodynamic principles [8, 12].

3.1. Pharmacokinetics: Simulating Drug Concentration Over Time

A successful therapy depends not only on a drug's efficacy but also on its concentration at the site of infection over time. This dynamic is known as pharmacokinetics (PK). Instead of assuming a constant drug level, our model simulates a concentration curve, $D(t)$, which rises after application and then gradually decays. This behavior is governed by two parameters:

• Time to Maximum Concentration ($T_{max​}$): This parameter defines the time in days for the drug to reach its peak concentration within the plant's system after administration ^[8].
• Drug Half-life ($t_{1/2}$​): This represents the time required for the drug's concentration to decrease by half. We defined distinct half-lives for CTX and oxytetracycline ($t_{1/2}$ = 100 and 200 days, respectively) to reflect the prolonged stability required for treatments in perennial plants. The value for oxytetracycline was chosen to align with the long-lasting effects of commercial trunk injection products like ArborBiotic™, which are applied annually [7,8].

3.2. Pharmacodynamics: Modeling the Dose-Response Relationship

Once the systemic drug concentration is known, we must model its effect on the CLas bacteria. This relationship, known as pharmacodynamics (PD), is rarely linear. We implemented a standard sigmoidal dose-response curve using a Hill function, which captures the threshold-like behavior of most drug interactions. This function is defined by:

• Half-effective Concentration ($EC_{50}$): This is the most critical PD parameter, representing the drug concentration needed to achieve 50% of its maximum effect. It defines the potency of the drug. In our model, we explored the hypothesis that a specifically designed antimicrobial peptide, such as CTX, could achieve a higher potency than a broad-spectrum antibiotic like oxytetracycline. To test this scenario, we set a plausibly lower $EC_{50}$ for CTX, based on potencies reported in the literature for other targeted antimicrobial peptides.
• Hill Coefficient (n): This parameter determines the steepness of the dose-response curve. A value of n=2.0 was chosen to model a sharp, switch-like activation, meaning the therapeutic effect engages strongly once the drug concentration approaches the $EC_{50}$ threshold.

The model distinguishes between the mechanisms of the two drugs. The efficacy of CTX is represented by a bactericidal term, Killing Efficacy (E_kill), which directly increases bacterial mortality. In contrast, oxytetracycline's effect is modeled as Growth Inhibition (E_inhib), a bacteriostatic term that reduces the bacterial replication rate. The maximum possible effect for either drug is modulated by the killScale parameter, allowing us to adjust the overall therapeutic strength.

3.3. Oxytetracycline-Specific Dynamics: Inhibition and Clearance

For oxytetracycline, we introduced an additional dynamic to better simulate its bacteriostatic mechanism. While the primary effect is to halt bacterial growth ($E_{inhib}​$), a static population would still occupy the phloem. Therefore, we included an Active Clearing Rate ($E_{clear}$​). This term models the host plant's ability to naturally remove weakened or non-replicating bacteria ^[10]. This clearing rate is directly proportional to the inhibition effect ($E_{\text{clear}} = \alpha_{\text{clear}} \cdot E_{\text{inhib}}$.), ensuring that as the antibiotic becomes more effective at stopping growth, the system also becomes more efficient at clearing the pathogen ^[7].

Parameter Summary

4. Fundamental Constants and Assumptions

While many parameters in our model are designed to be adjustable to reflect biological variability, a set of fundamental constants provides the stable foundation upon which the simulation is built. These fixed values, defined in the constants.h file of our source code, establish the "laws of physics" for our virtual environment. They represent core assumptions about the system's capacity, numerical stability and computational framework, ensuring consistency and reproducibility across all simulation runs.

4.1. Infection Initiation and Spread

To simulate the progression of the disease, we first needed to define how an infection is established and how it moves from one cell to another.

• Initial Infection Load ($I_0$​): This constant sets the starting bacterial concentration (set to 0.1) for the single cell where the simulation begins. This "patient zero" approach allows us to observe the emergence of spatial patterns from a controlled origin point, mimicking the initial infection by an insect vector.
• Spread Infection Load: This parameter defines the fixed amount of bacterial load (0.05) that is transferred from an infected cell to an adjacent healthy cell during a propagation event. By keeping this value constant and relatively low, we model the slow, cell-to-cell spread characteristic of CLas, where the pathogen must overcome physical barriers to advance through the phloem.

4.2. Host Defense Activation Threshold

As discussed in the biological dynamics section, the plant's callose response is not linear but rather a threshold based mechanism. The activation of this defense is governed by two fixed constants that define the underlying Hill function:

• Callose Signal EC50​ ($EC_{50}^{C}$​): This is the activation threshold for callose production. The value is set to 0.5, meaning that the defense response is triggered at 50$\%$ of its maximum rate only when the local bacterial load reaches half of the cell's carrying capacity. This ensures the plant does not overreact to very low pathogen levels but mounts a robust defense once a significant infection is established.
• Callose Hill Coefficient ($n_C$​): A coefficient of 2.0 was chosen to create a steep, switch-like response curve. This models a decisive activation of the defense mechanism once the bacterial signal crosses the EC50​ threshold, consistent with the principles of PAMP-triggered immunity (PTI), where a sharp response is mounted upon pathogen recognition.

4.3. Therapeutic and Numerical Stability Assumptions

Finally, certain constants were established to ensure specific model behaviors and computational robustness.

• Oxytetracycline Active Clearing Coefficient ($\alpha_{clear}$​): As described in the therapeutic module, our model for oxytetracycline includes an "active clearing" effect, where the host helps remove bacteria whose growth is inhibited. This constant defines the fixed efficiency factor (0.015) for this process $\alpha_{\text{clear}}$ and $E_{\text{clear}} = \alpha_{\text{clear}} \cdot E_{\text{inhib}}$. This value was chosen to be small, representing a slow, secondary process compared to the main infection dynamics and it represents a core assumption that there is a slow but steady host-mediated removal of the static pathogen under bacteriostatic treatment.
• Numerical Extinction Threshold: Set to $10^{−3}$, this value represents a computational cutoff. Any bacterial concentration within a cell that falls below this threshold is considered zero. This prevents the model from spending computational resources on biologically insignificant pathogen levels and avoids potential numerical floating-point errors, ensuring simulation stability.

Parameter Summary

The following sections present the findings generated by our computational model. We begin by conducting a validation, comparing the model's output against experimental data to establish its biological plausibility. With the model validated, we then leverage it as a tool to compare the efficacy of CTX peptide against conventional treatments. Finally, we present a sensitivity analysis to identify the most influential parameters governing the disease's behavior.

1. Model Validation: Reproducing Cellular Mechanisms

Before using our model to simulate drugs, it was essential to validate its behavior against real world biological data to ensure that our "virtual laboratory" captured the key characteristics of Greening. After defining the baseline parameter set for our control simulations , we compared our model results with experimental results at three distinct scales: spatial (cellular mechanisms), temporal (long-term progression), and dynamic (emerging patterns).

1.1 Simulation Baseline Parameters

The validation results presented in this section were generated using a single, consistent set of baseline parameters. These values, chosen from the ranges defined in the "Model Parametrization" section, represent a biologically plausible scenario for the progression of Citrus Greening. The specific values used for these simulations are summarized in the table below. The fundamental constants of the model remained fixed as described previously.

1.2. Spatial Validation: Reproducing Cellular Mechanisms

Our model was designed as a 'virtual laboratory' to translate complex biological interactions into computational rules. This section demonstrates how our model successfully captures the fundamental dynamics of Greening, starting with the spatial "tug of war" between the pathogen CLas and the plant's defenses at the cellular level.

Our first objective was to verify if the model could reproduce a key, counter-intuitive observation described in the literature: the spatial displacement of the callose response. The figure below places the biological evidence that inspired our model side-by-side with the emergent spatial patterns generated by our simulation.

Our model is grounded in a key observation from Bernardini et al. (2022) [1]: phloem pores in CLas-infected cells are, counter-intuitively, significantly more open than in uninfected cells (Fig. 1A, 1B). This suggests the bacterium actively manipulates the host's primary defense: callose deposition. Our simulation was built to test this hypothesis.

The results strongly corroborate these experimental findings. The quantitative analysis of our simulation (Fig. 2D) shows that "Infected" cells have a drastically lower mean callose concentration than neighboring "Uninfected" cells. This aligns with the experimental findings at both the physical level (open pores) and the genetic level, where. found that genes for callose synthesis were downregulated in infected tissues [1]. According to the permeability mechanism implemented in our model ($P(t) = e^{-\lambda_C \cdot C(t)}$, Fig. 2B), a lower callose concentration ($C$) directly results in higher permeability (P), that is a more open pore. Therefore, our model not only reproduces the experimental observation but also provides a clear mechanistic explanation for it, serving as a primary validation of our model's core logic.

Furthermore, the simulation reveals the spatial dynamics of this interaction, the "dynamic interplay" in action. The heat maps (Fig. 2A) visualize a front of infection (red) being contained by a "wall" of callose (yellow) that forms in adjacent sites. This is confirmed by the radial profile (Fig. 2C), which shows the callose peak occurring at a distance from the infection's center, forming a "ring of defense" around the compromised area. This result demonstrates that our model successfully captures the sophisticated spatial displacement of the immune response, a key feature of the CLas-citrus interaction [1].

1.3. Temporal Validation: Simulating Long-Term Disease Progression

After establishing that our model reproduces local cellular mechanics, the final test of its validity was to verify if these rules, when simulated over a long period, would generate a disease progression curve that matches field data. In other words, are the local interactions sufficient to explain the macroscopic disease dynamics?

Comparative Analysis of Disease Phases

Figure 3 shows a remarkable agreement between our model's curve and the normalized field data. The analysis reveals that our model successfully captures the three key phases of the disease:

1. Latency Phase (0 to ~4 months): The model exhibits an initial period of very slow bacterial growth, with the infection remaining at low levels. This perfectly corresponds to the latency period observed in the field, where the plant does not yet show visible symptoms.
2. Exponential Growth Phase (4 to 10 months): The simulation shows a sharp increase in infection, which aligns with the temporal window where disease symptoms begin to manifest (gray area). The data from Coletta-Filho et al. (2010)[2] fall directly on our curve's upward trajectory, validating the timing of the disease progression.
3. Saturation Phase (after ~1 year): Finally, answering our validation question, the simulated infection curve reaches a stable plateau at $\bar{I}(t) \approx 0.4$. This equilibrium state aligns perfectly with our "anchor point" extracted from Vasconcelos et al. (2020)[3], which represents the chronic infection level in the field. This agreement demonstrates that the model successfully ilustrates the disease's final saturation state.

Validation Methodology and Data Normalization

To validate our model's temporal trajectory, we selected two key studies that characterize Greening progression at different scales: the acute phase (months) and the chronic phase (years).

• Acute Phase Data (Coletta-Filho et al., 2010)[2]: This seminal study characterized the acute infection dynamics in young orange trees. In the experiment, healthy plants were infected via grafting, and the CLas bacterial load was quantified at regular intervals (30 to 240 days) using qPCR. Their results revealed a clear phased progression: a latency period with a low, constant bacterial population for up to 90 days, followed by a rapid growth phase where the concentration increased over 10,000-fold, peaking around 180-240 days. We used this data because it provides a precise temporal "ruler" for disease establishment, making it the ideal benchmark to validate the timing of our model's latency and exponential growth phases.
• Chronic Phase Data (Vasconcelos et al., 2020)[3]: To validate the long-term behavior, we used the study from Vasconcelos and collaborators (2020) as our "anchor point". This work quantified the total CLas population and its distribution in 2.5-year-old 'Valencia' orange trees with chronic infection. They found that the bacterial load stabilizes at approximately $8.3 \times 10^6$ to $9.9 \times 10^6$ cells/g in canopy tissues (leaves and branches). We used this data because it represents the stable equilibrium state of the disease in the field.

To rigorously validate our model's temporal trajectory, we adopted a two-step validation approach. First, we used the data on the chronic state of the disease (Vasconcelos et al., 2020)[3] as a fundamental biological constraint to calibrate our parameter space. We ran multiple simulations and selected only the parameter sets that could achieve the correct infection plateau. Second, and most critically, we then tested whether this calibrated parameter set could reproduce, without further tuning, the data from an entirely independent study on the acute infection phase (Coletta-Filho et al., 2010)[2]. The fact that the parameters which explain the end of the disease also described the beginning constitutes a strong cross-validation, demonstrating that our model captures the fundamental dynamics of Greening progression.

1.4. Dynamic Validation: Emergence of Spatial Patterns

While the temporal analysis (Figure 3) validates the macroscopic progression of the disease, the beauty of a cell-based model lies in its ability to simulate spatial dynamics. Figure 3 presents an animation of our control scenario simulation, visualizing the "cellular battle" between infection and defense in action over time.

The animation reveals that the infection does not spread as a uniform battlefront; instead, it advances slowly during the initial phases, forming "branches" that explore paths through the tissue. In response, the callose defense forms a "wall" that constantly attempts to contain the pathogen's advance.

As time progresses, these initial infection clusters expand and merge until the pathogen dominates most of the tissue, leaving behind only "islands" or isolated clusters of healthy cells, as seen in the late stage of the simulation. This final heterogeneous pattern, with pockets of healthy tissue remaining, is a key qualitative validation of our model. It visually represents the biological reality that CLas is not uniformly distributed throughout the plant's phloem [1], confirming that the local rules we programmed are sufficient to generate the complex infection patterns observed in nature.

1.5. Conclusion of Model Validation

This multi-level validation is a fundamental step. It demonstrates that the local, short-term interactions we programmed are sufficient to explain both the macroscopic, long-term dynamics and the complex spatial patterns of the disease. This gives us the confidence to, in the next section, use our model as a tool to test our therapeutic interventions.

2. Analyses of Therapeutic Interventions

Having established the validity of our model against experimental data, we then leveraged it as a tool to compare the efficacy of our CTX peptide against a conventional antibiotic treatment and a no-treatment scenario. The simulations were initiated from a chronic infection state and run for an extended period to observe the long-term outcomes of each intervention.

2.1. Therapeutic Parameters and PK Profiles

To simulate these treatments, we first defined their distinct pharmacokinetic (PK) and pharmacodynamic (PD) properties. The specific parameters governing the drugs' behavior and potency are detailed in Table 3. The PK profiles in Figure 4 show how the concentration of each drug evolves over time, a critical factor that dictates its therapeutic effect[5].

2.2. Comparative Efficacy Analysis

With these distinct mechanisms defined, we applied each treatment to the chronically infected system. The comparative results, shown in Figure 6, clearly demonstrate the superior efficacy of our bactericidal approach.

Analysis of Therapeutic Outcomes

The direct comparison of the three simulated scenarios, visualized in Figure 6, reveals the distinct long-term impact of each therapeutic strategy on both the pathogen and the host's immune response.

• Control Scenario: The top panel serves as our baseline. As validated previously, the infection progresses to a high, chronic saturation point ($\bar{I}(t) \approx 0.4$). This sustained bacterial presence triggers a continuous, high-level callose response ($\bar{C}(t)$), representing the chronic phloem plugging that characterizes the disease's symptoms.
• Oxytetracycline (Bacteriostatic Action): The middle panel demonstrates the effect of a conventional bacteriostatic antibiotic. The application of oxytetracycline causes a sharp drop in the infection level, and correspondingly, the callose level also decreases, suggesting a temporary alleviation of phloem-plugging symptoms. However, due to its bacteriostatic nature, the treatment only suppresses the bacterial population. Once the drug's concentration falls below a therapeutic threshold, the infection begins to rebound. As the pathogen population regrows, the host's defensive callose response is re-triggered, and the callose level rises back towards its chronic, damaging state. Interestingly, in this rebound scenario, the infection reaches a new peak ($\bar{I}(t) \approx 0.8$) that is significantly higher than its initial chronic state, indicating a potential dysregulation of the host response.
• CTX Peptide (Bactericidal Action): In stark contrast, the bottom panel shows the definitive impact of the CTX peptide. Its bactericidal action causes a rapid and profound decline in the bacterial population. As the infection is eliminated, the trigger for the host's immune response vanishes. As a result, callose deposition ceases completely, with its level ($\bar{C}(t)$) returning to zero. This represents not just the removal of the pathogen, but the full recovery of the phloem from the diseased state.

2.3. Conclusion: A Curative vs. Suppressive Strategy

The comparison between the treatments illustrates the fundamental difference between a suppressive and a curative strategy. To explore the potential magnitude of this difference, we ran a specific "what-if" scenario based on plausible parameters for a targeted peptide. Within this simulated scenario, a curative outcome was achieved with a CTX dose ten times lower than that of oxytetracycline. While not a quantitative prediction of the real-world value, this result provided us with a justification interesting enough to guide our experimental design (see the full results on our Engineering page), justifying the decision to begin our in-planta tests with a correspondingly low dose.

To complement the temporal analysis, we generated animations of the spatial grid for each of the three scenarios. These visualizations provide an intuitive understanding of how each treatment impacts the disease at the tissue level, translating the graph curves into a concrete spatial outcome. Yellow represents callose, red infection and black a healthy site.

The animations in Figure 6 provide a powerful visual confirmation of the results seen in the temporal graphs. The Control simulation depicts the unchecked invasion, culminating in a grid dominated by chronic infection. The Oxytracycline simulation is particularly revealing, showing an initial regression of the infection followed by a clear rebound, where isolated pockets of surviving bacteria re-colonize the grid. In stark contrast, the animation for the CTX peptide shows a rapid and definitive clearance of all pathogen clusters until the grid is completely healthy. This visualization provides the spatial proof of complete eradication, representing a curative outcome.

3. Sensibility of parameters

To understand which biological factors have the greatest impact on disease progression and to test the robustness of our 'virtual laboratory,' we performed a sensitivity analysis. In this process, we systematically varied key parameters governing the model's behavior and observed the effect on the temporal and spatial outcomes of the simulation.

We first analyzed the parameters that control the kinetic "arms race" between the pathogen and the host: the bacterial growth rate ($r$), the callose production rate ($\alpha_C$), and the bacterial propagation rate ($\beta$).

The analysis reveals how each parameter governs a different aspect of the disease:

• Bacterial Growth Rate ($r$): Panels (a) and (b) show that the model is highly sensitive to the intrinsic growth rate. Higher $r$ values lead to a much more severe chronic state (higher saturation level), suggesting that the aggressiveness of replication is a primary factor dictating the final severity of the disease.
• Callose Production Rate ($\alpha_C$): Panels (c) and (d) illustrate the "tug-of-war" between pathogen and host. A strong callose response (high $\alpha_C$) is able to contain the infection at a low chronic level. In contrast, a weak defense allows the infection to completely saturate the tissue, validating callose as a crucial containment mechanism.
• Bacterial Propagation Rate ($\beta$): Panels (e) and (f) show that the $\beta$ parameter primarily controls the speed at which the disease spreads. Higher values lead to a much faster saturation, indicating that $\beta$ determines how quickly the plant reaches its chronic state.

Next, we investigated the impact of the Defense Signaling Radius ($S_r$), which defines the host's spatial defense strategy.

The results show this parameter is decisive in shaping the disease's spatial footprint:

• A small signaling radius ($S_r = 1, 3$) leads to a reactive and ineffective defense, allowing the pathogen to spread in a diffuse, uncontrolled manner.
• In contrast, a large signaling radius ($S_r = 9, 12$) allows the plant to mount an extremely effective "quarantine," strangling the infection within very small, isolated clusters.
• We observed that an intermediate signaling radius ($S_r = 6$), our baseline value, generates a complex, heterogeneous pattern that is qualitatively consistent with the sparse infection distributions seen in nature.

In summary, our sensitivity analysis demonstrates that the model behaves in a biologically robust manner. It reveals that the final severity of the disease is primarily determined by the "arms race" between bacterial replication ($r$) and the strength of the host defense ($\alpha_C$), while the speed of progression is governed by $\beta$. Finally, the spatial pattern of the infection is dictated by the host's signaling strategy, $S_r$. This understanding reinforces our model's validity and provides insights into which factors are most critical in the CLas-citrus interaction.

Our computational model was not an add-on, but a powerful tool that guided our project from initial strategy to experimental design, allowing us to make data-driven decisions and optimize our efforts.

1. Defining the Strategy to Combat Greening

The starting point of our project was a strategic question: what is the most effective approach against CLas? A bacteriostatic suppression or a bactericidal eradication? To answer this, our model served as a virtual laboratory, testing both scenarios. The results were unequivocal: the simulation demonstrated that the bacteriostatic approach led to a severe disease relapse, making it an unsustainable strategy. In contrast, the bactericidal approach, like that of our CTX peptide, proved to be curative, eradicating the pathogen and allowing for the full recovery of the tissue. This conclusion was a cornerstone of our project.

2. Validating the Biological Plausibility of Our System

To ensure our conclusions were reliable, we first made sure the model was biologically plausible. The model's first contribution was to successfully simulate and reproduce the "spatially displaced callose response," a complex mechanism of the pathogen-host interaction. This validation gave us the confidence that our "virtual laboratory" was a sufficiently faithful representation of reality to guide our decisions.

The model's most direct contribution was bridging the gap between theory and practice in the laboratory.

1. Guiding the Experimental Design: A central question was what dose to begin our in-planta assays with. Our model provided immense practical value here. By parameterizing it with plausible values from the literature, the simulation indicated that a curative effect could be achieved in the low-milligram range. This insight did not provide a final, optimized value, but it gave us the confidence to focus our experimental tests on a promising and manageable range, starting with 15 mg (see the full results on our Engineering page). The model thus served as a strategic tool to optimize our experimental design and resources.
2. Capabilities of the Modeling Tool: Our model was built as a robust platform to answer multiple questions. Using it, we can:
• Simulate various infection scenarios and plant physiological responses.
• Explore different dosage scenarios to help identify promising ranges for reducing the spread of CLas.
• Generate infection maps to identify the most vulnerable phloem regions.
• Estimate infection percentages and visualize emergent patterns, informing field intervention strategies.

In summary, our model validated our strategy, gave us confidence in our biological understanding, and most importantly, provided strategic guidance that made our experimental work more focused, efficient, and with a higher probability of success.

We recognize that every model is a simplification of the complex biological reality, designed to help us understand fundamental patterns in nature. The honest identification of our work's limitations is not only good scientific practice but also the map that guides our next steps. We present here the main limitations of our current model and the future vision they inspire, establishing a clear path for our work and that of future teams to continue to advance.

1. Capturing Systemic Dynamics

• Current Limitation: Our main simplification was to model the phloem as a 2D grid. This was a deliberate choice to isolate and understand local cellular dynamics. However, this approach does not capture the long-distance vascular transport network.
• Future Work: A possibility could be to adapt this into a graph.

2. Feeding the Model with More Experimental Data

• Current Limitation: Many of our parameters were estimated based on literature data to ensure a qualitatively correct behavior. Our validation confirmed the model's plausibility, but not its quantitative predictive power.
• Future Work: We intend to pursue two fronts to increase the model's rigor:
• Experimental Collaboration: Establishing partnerships with laboratories to obtain empirical measurements of key parameters, such as the in planta replication rate of CLas.
• Predictive Validation: Implementing a rigorous validation where the model would be calibrated on one dataset and then tested against a completely independent dataset to verify its predictive capability.

3. Advancing the Therapeutic Module: Spatially-Resolved Pharmacokinetics

• Current Limitation: Our model assumes a uniform systemic concentration of our CTX peptide after application.
• Future Work: We plan to develop a spatially-resolved pharmacokinetic/pharmacodynamic (PK/PD) module. This new module would simulate how the peptide diffuses and is transported through the xylem and phloem, allowing us to optimize not only the dose but also the method and site of application for maximum efficacy.

4. The Legacy of Our Project

Beyond improving the current model, our work lays the foundation for two exciting future directions:

The modular architecture of our model was designed to be adaptable. The core logic that simulates the interaction between pathogen, host defense, and therapy can be parameterized for other diseases. Our work can serve as a foundational model for the scientific community, allowing for the development of "virtual laboratories" for other important diseases that follow similar dynamics.

Our model not only answered questions but also generated new ones. The prediction that bacteriostatic treatments may lead to a more severe disease relapse is a fascinating and testable hypothesis.

• Future Work: This prediction directly inspires a new line of experimental research in our laboratory. The next cycle of our project could focus on validating or refuting this hypothesis, demonstrating the power of the complete engineering cycle: the model informs the experiment, and the experiment refines the model.

Our project is built on the principles of open science and reproducibility. On this page, you will find all the necessary resources to explore, reproduce, and build upon our computational work, including the full source code for our simulation and additional supplementary materials that provide deeper discussions about our model and some curiosities of this type of science.

Code Repository

All the source code for our model is publicly available on our team's GitHub repository. The repository contains the complete simulation engine written in C++, along with the Python scripts used for data analysis, plotting, and video generation presented on our Results page.

If you don’t have familiarity with C++ we translate our model for Python, but it’s important to understand that the process in C++ is faster than Python, and because of that we program in C++.

We invite you to explore, download, and run our code.

Code Repository

Are you Curious? Dive Deeper!