Rice blast occurs in major rice regions of the world, with the more severe rice areas concentrated in Asia and Africa. Rice blast can occur in all growth stages and parts of rice, the most serious of which are leaf blast and neck blast. In disease-endemic areas or diseased fields, about 10%~30% of rice yield is reduced or even extinct every year [1]

Major grasses (maize, wheat and rice) together account for about half of global caloric intake. Unfortunately, our global dependence on grasses has inadvertently promoted some of the most destructive fungal plant pathogens, such as Magnaporthe Oryza sativa, the pathogen of rice blast.

In order to resist the damage of viruses to plants, we designed saRNA-based plant vaccines, and modeled the processes of virus infection, vaccine transport and delivery, efficacy time, antigen-antibody matching, etc. The whole modeling adopts a cross-scale model, from the macroscopic infection process to the microscopic molecular transport, combined with the neural network model and finite element method, and compared with the experimental results, to provide explanations and feasibility demonstrations for biological phenomena. The entire model is divided into four parts:

1.1 Wind Field Large Eddy Simulation (LES)

Meteorology is a huge factor affecting the whole process of plant growth, and there are various air currents around the world, as shown in the figure below [2]. Since the transmission mode of rice blast spores is mainly environmental transmission (such as airflow, etc.), it is necessary to simulate the airflow field eddy in the planting environment. The Yuquan Campus of Zhejiang University (longitude: 120.123077, latitude: 30.263842) where the ZJU-China team is located was used as an example to simulate the rice planting environment.

Because it is large-scale meteorological flow, then, large eddy simulation is used to directly calculate large-scale fluid pulsation, and small-scale vortex is filtered. For the fundamental Navier-Stokes tensor equation:

where $u_i$ represents the magnitude of velocity, $x_i$ represents the spatial position of the fluid, $ f_i $ represents the mass force, $\rho$ represents the density of the fluid, $p$ represents the pressure, $t$ represents time, and $\nu$ represents Kinematic viscosity.

The method of small-scale vortex filtration is to perform low-pass filter of the Navier-Stokes equation, which decomposes the fluid velocity into large-scale motion velocity (solvable scale pulsation) and small-scale pulsation velocity (sublattice-scale pulsation):

where $u(\boldsymbol{x},t)$ represents the actual velocity, $\tilde{u}(\boldsymbol{x},t)$ represents the low-pass pulsation, and $u''(\boldsymbol{x},t)$ represents the residual pulsation. The integral filtering equation can be determined by the following expression:

where $G(\boldsymbol{x} - \boldsymbol{\xi})$ is the filter function, and the integral filtering of the cube of length $l$ is performed, then the one-dimensional box integral filtering function can be expressed as:

In addition to box filtering, there are also gaussian filters and other forms. Considering the uniform filtering process, in the physical space, the filtering operation and the partial conduction operation can be interchanged. Substitute the above form into the Navier-Stokes equation, we can derive the filtered Navier-Stokes equation:

where $\widetilde{(\cdot)}$ represents filtered variable, and $\tau_{ij} = \tilde{u}_i\tilde{u}_j - \widetilde{u_i u_j}$ is the sublattice stress. Sublattice stress is the momentum transport between filtered out small-scale pulsations and solvable-scale turbulence. To achieve large eddy simulation, it is necessary to construct a closure of sublattice stresses.

Assuming that the small-scale pulsations filtered out by the isotropic filter are locally equilibrium, that is, the energy transfer from the unsolvable scale pulsations of the solvable scale term is equal to the turbulent energy dissipation, the sublattice Reynolds stress mode in the vortex-viscous form is [3]:

where $\nu_t$ is the sublattice vortex viscosity coefficient, $S_{ij}$ is the strain rate tensor, and the figure below is the large eddy simulation results of the environment of Yuquan Campus of Zhejiang University in a defined physical space.

Assuming that the wind flows from the cold source to the hot environment, the West Lake is considered as the cold source, so the setting condition is that the air flows from east to west. From the simulation results, it can be seen that in the part where the mountain and architecture areas are not encountered, the air flow has relatively small vortex generation, while in the part with buildings and mountains, the air flow has a relatively large and unstable vortex state, and the streamline line is also deflected, so the planting environment here has the characteristics of turbulent atmosphere.

1.2 Turbulent diffusion model of plant viruses

Consider a large, homogeneous field, planted to the rice which is genetically uniform, and a single infection of the plant pathogen is introduced to the field. Consider further that the theory of turbulent diffusion governs the dispersal of windborne spores of most foliar pant pathogens which suggests that spread becomes more efficient as the area with diseased plants expands. This should result in epidemics with increasing frontal velocity (dispersive epidemic waves).

Under ideal circumstances, the direction from which the wind moves through the rice is random. Then, the fact leads to a focal epidemic whose front will expand radially in a wavelike manner from the site of initial infection [4].

Let $U(r,t)$ be the diseased rice tissue (e.g., the number of diseased leaflets per plant) at radial distance $r$ and time $t$. Based on the simulation of the airflow above (see Section 1.1 for details) with the theory of Richardson for turbulent atmospheric diffusion predicts that the diffusion coefficient is $K(r) \propto r^{\frac{4}{3}}$. And considering the disease intensifies follows the Vanderplank logistic equation with apparent infection rate $\gamma$, then the differential equation is:

where $K(r) = C_1r^{\frac{4}{3}}$, $C_1$ is a constant, then:

However, the equation (10) cannot be solved analytically. If we only focus on the front of the epidemic, the differential equation can be solved approximately by changing the logistic term $\gamma U (1-U)$ into linear term $\gamma U$:

solving the equation, together with the boundary condition $U(r,0) = \delta(r)$ and $U(\infty, t) \rightarrow 0$:

The table above shows the changes in the proportion of infections measured with the time series at a defined place $r$. According to the quasi-steady-state frontal approximation, the curves are fitted using the sets of data, and the Trust-Region Dogleg method is used to calculate the nonlinear system of equations with random initial iterations, and the constant values and iterative data are obtained in the following table (only convergent solutions are shown):

It can be found that the convergent solutions are all the same convergent solution in the error range, other initial iteration points do not make the results converge or conform to practical requirements, and the larger the value of the random point (the more it deviates from the initial stability domain), the less it converges.

From the left graph, it can be seen that the data can match the curve only if it meets the quasi-static short time approximation, and the curve deviates after exceeding the quasi-static short time.

1.3 Plant virus pattern dynamics

In order to construct a mathematical model consisting of activator and inhibitor, assuming $V(r,t)$ as the dead or recovered rice tissue (i.e., that tissue cannot infect other plants) at radial distance $r$ and time $t$ , the equation is:

where $K_2$ is the diffusion coefficient of the infection-free rice, $\tilde{G}(U,V)$ is the function controls interactions between nfected and uninfected plants. Considering a small area of the field, the diffusion coefficient can be regarded as a constant parameter $K_1$, then the equation is:

Combine equations (13)(14) and take a polynomial function approximation to the function $\tilde{F}$ and $\tilde{G}$, assuming that the uninfected tissue will get infected by contacting with the other two infected tissue, then the form is consistent with pattern dynamics.

The whole process of the spread of M. Oryza sativa is shown in the figures above, the patterns start from small random position, and then spread to the whole area, and the speed of diffusion first increases and then decreases.

1.4 Plant virus infestation process

The rice blast is mainly in the form of conidia and hyphae under certain temperature and humidity conditions relying on airflow, through airborne transmission to plant leaves and other aerial parts, in the presence of water film and suitable environmental conditions, the viscosity of the conidia tip makes it adhered to the surface of aboveground tissues such as rice leaves. Then, under the premise of not needing external nutrients, the conidia only germinate through the endogenous nutrients of the spores to form bud tubes, and then the bud tubes specifically differentiate into attached cells with melanin. Then the attachment cells of melanin form infection plugs and penetrate the stratum corneum and superficial cell wall of some plant tissues above ground, invading epidermal cells or elongated cells. After invasion, rice blast bacteria grow in rice cells, and the hyphae rapidly expand and spread, and then infect adjacent epidermal cells and penetrate deep into mesophyll cells.

After the infection is completed, hyphae and conidia are produced in the infected cells, and are directly released into the air from the stomata or penetrating the plant epidermis, and new conidia are produced at the apex of the conidia peduncle, and mature conidia are released into the air, and again spread to the leaves of other rice by means of airflow propagation to cause secondary infection [5].

2.1 Physical Shapes and Entropic Elasticity of RNA in Fluids

Consider random walk of saRNA to describe its physical morphology in rice vascular. Assuming the polymer chains can be decentralized into plenty of rigid segments, denoted Kuhn length $a$. It is assumed that the physical morphology of saRNA is determined by end-to-end distance $R$, i.e., the distance that the end of the RNA segments leaves the initial end after random walk $N$ times. Using the averaging method, the expected value of $R$ is:

where $x_i = \pm a$ is the displacement of step $i$, then considering the variance of random walk:

the expansion of equation (16) is:

Due to the randomness of the walk, the walk at step $j$ is not correlated with the $i$ step, then the second term of equation (17) is 0, and the variance is determined entirely by the autocorrelation displacement:

according to the results of one-dimensional random walk [6], it can be generalized to three-dimensional form:

the parameters $\mathcal{N}$ and $\kappa$ are determined by normalization and variance:

the random walk chain distance distribution consisting of Kuhn length with $N$ step length $a$ is obtained:

From the results, it can be seen that ideally, the probability of conformation is the largest when $R$ is small (the probability of center position is the highest in the figure), that is, RNA molecular chains tend to form more uniform clusters, while fully elongated long chains are not easy to appear, and the simulation schematic diagram of the two cases based on molecules is as follows:

Appendix

The theory that RNA segments can be discretized into plenty of rigid segments is based on the theory of persistence length of RNA, assuming that within the persistence length, the molecular segments basically take the form of straight segments. The persistence length is the characteristic scale of the tangent vector correlation function attenuating along the chain segment:

where $\boldsymbol{t}(x)$ represents the tangent vector, $x$ represents the position along the segment in a chain segment of length $L$, and $\xi_p$ is the resident length, it should be noted that the exponential decay rate only applies to sufficiently small spacings $\left| s-u \right|$. By definition, the end-to-end distance is:

then the end-to-end distance variance is:

It can be seen that the relationship between the Kuhn length $a$ and the resident length $\xi _p$ is $a = 2 \xi _ p$.

The physical morphology of RNA is not only affected by thermodynamics, but also by the electrostatic action of the spatial structure, as well as the electrical action of ions around the liquid during vascular transport (satisfying the Poisson-Boltzmann equation), considering the above influencing factors, the ball and stick model of a physical form of RNA is shown below [7]

2.2 Wiener process in vascular bundle

Based on the investigation of its physical form in Section 2.1, if the saRNA molecule is regarded as a whole, it meets the requirements of the Wiener process for molecular scale of 0.1~10 microns. Consider RNA molecule with coordinates $\boldsymbol{q}$ and momentum $\boldsymbol{p}$ suspended in a liquid composed of solvent particles:

where the $V_{FP}$ is the interaction of solvent particles and particles, the $V_{FF}$ is the interaction of solvent particles and solvent particles, and the $V_{ext}$ is the external field potential. Leaving aside the macroscopic flow of the fluid, assuming that there is a harmonic oscillator potential $V_{FP}(\boldsymbol{q} - \boldsymbol{q}_i) = \frac{1}{2}m_i\omega_i^2(\boldsymbol{q} - \boldsymbol{q}_i)^2$, and assuming that the interaction potential between solvent molecules is zero potential energy, and all solvent molecular masses are the same $m$, and $M \gg m$, then the Hamiltonian can be simplified to: $ H(\boldsymbol{p}, \boldsymbol{q}, \{\boldsymbol{p}_i, \boldsymbol{q}_i\})$ $=$ $\frac{\boldsymbol{p}^2}{2M} $ $+$ $\sum_{i}\frac{\boldsymbol{p}^2}{2m_i}$ $+$ $\frac{1}{2}\sum_i \omega_i^2(\boldsymbol{q} - \boldsymbol{q}_i)^2$ $+$ $V(\boldsymbol{q}) $

according to the momentum theorem and Hamiltonian equations:

From the Duhamel integral, it can be seen that the trajectory of the solvent molecule is:

Integral by parts to equation (27):

Then the Hamiltonian equation for RNA clusters is:

where $ (\boldsymbol{q} - \boldsymbol{q}_i) $ can be obtained by transposition from equation (28):

let:

Substitute the above equation into (29)(30), we can derive the equilibrium dynamical equation:

It can be seen from the above kinetic equation that without considering the interaction between solvents, $K(t)$ reflects that the solvent molecules only play a buffering role, remembering the effect of RNA particles on the solution at the past moment and delaying this effect to the RNA particles; $\xi(t)$ is related to the initial states of the system, which are unknown and can be treated as random variables. Assuming that the solution is in equilibrium, $\xi(t)$ conforms to the Gaussian distribution, then there are:

Taking $K(t) = 2 \mu \delta(t)$ gives the Langevin equation for RNA, where $\xi(t)$ is white Gaussian noise. Since the movement of RNA is a superposition of random motion and transpiration tension, if the influence of one of these factors is studied, it is necessary to characterize it with the assistance of another equivalent tracer particle, because the expression of fluorescent protein is related to RNA, and to a certain extent, the two are similar in scale, so fluorescent protein can be used as a "conjugated particle" to approximate the influence of RNA random motion factors. For solving the above stochastic differential equation, the following assumptions are met according to the experiment:

Since the mass of the fluorophore is small, the kinetic equation after ignoring the inertia term can be reduced to:

Solving the above equation, we can derive $q(t) = W_t$, i.e., the solution to the differential equation for the stochastic process is the Wiener process. For a stochastic $X_t$ process, if it has a Wiener process with an offset of $f(X_t, t)$ and a standard deviation $\sqrt{2D(X_t, t)}$, the corresponding Langevin equation is:

For the random process to take the value at the left end of the interval, the Ito- Langevin equation for fluorescence motion can be obtained by using the Ito integral:

The trajectory can be obtained by computer simulation of the above process, assuming that the fluorescence starts from the black point to the red point terminate, and the trajectory color changes from the initial position white to the dark blue at the end point, which is consistent with the fluorescence diffusion of the experiment under statistical average, as shown in the figure below:

2.3 Transpiration tension in vascularity

Rice transpiration can be modeled according to energy conservation and mass conservation, which links transpiration to leaf-to-air vapor pressure difference:

where $\lambda$ is the latent heat of vaporization, $ET$ is the evapotranspiration rate, $R_n$ is the net radiation intercepted by the crop, $H_c$ is the sensible heat exchanged between the canopy and the air, $\rho$ is the density of the air, $C_p$ is the specific heat of air, $\gamma$ is the psychrometric constant, $g_t$ is the total canopy conductance to water vapor transfer, $D_c$ is the canopy-to-air vapor pressure difference. Since some parameters are difficult to measure, the parameters need to be transformed or equivalented [8].

In the process of modeling plant conducting tissues, the transpiration of plants will cause the difference in water content in the ducts from roots to leaf fragments, resulting in inorganic salt concentration gradients, while the TLS segment of RNA can bind proteins to drive RNA movement in the ducts with a concentration gradient. Not only that, transpiration itself also causes water to flow from the roots to the leaves, and can also carry RNA throughout the body. Since the scale of the simulation is macroscopic, a cross-scale model needs to be established to simulate the average effect of the Wiener process in Section 2.2.

In the simulation, the chemical potential difference and osmotic pressure difference caused by the concentration gradient can be equivalent to the duct pipe pressure difference, and the overall equivalent structural pressure difference and the laminar flow of the pipe triggered by the initial flow rate can simulate the process of RNA passing through the plant conducting tissue to the whole body.

In the figure above, a small scale approximation part along the axis length of a plant vascular bundle is taken, and due to the small axis length of the selected segment, it can be approximately equivalent to a straight segment for a curved plant vascular. Considering the fluid pressure differential drive and Brownian motion caused by transpiration tension, the saRNA transport is relatively stable. And the actual plant vascular is with pipe curvature and there is a certain ring in the duct, ring and pipe bending will affect the flow state of fluid particles, making the actual flow more complex, as shown in the following two figures, vascular bending and the existence of ring lines will delay the flow speed, considering gravity, will further increase the resistance of the transpiration tension of the plant from the root to the leaf, so the flow state in the vascular tube is similar to laminar flow, the flow rate and streamline are relatively stable.

In the experimental results, taking the results of 720bp GFP mRNA and 3951bp RUBY mRNA as an example, three TLS and a control group were used in the experiment, and Actin was used as the control group in each collection, and the fluorescence intensity was measured with qPCR to reflect the mRNA concentration of different spatial positions of the plants on the fourth day.

3.1 Replication and degradation dynamics

Triggering the replication of saRNA requires the replicase to recognize a specific site, so the saRNA that is not bound to the replicase is regarded as an inactive state $I_s$ and the saRNA that binds to the replicase is regarded as an activated state $A_s$, then there are:

The state function $p$ is used to digitize the active and inactive states $\varphi(I_{\text{saRNA}}) = 0; \varphi(A_{\text{saRNA}}) = 1$, and assume that the degradation products of saRNA are $\eta$, and the degradation products of proteins are $\xi$, $\alpha,\ \beta,\ \gamma,\ \omega,\ c_1,\ c_2$ are constants:

A ternary function is used to describe a random system, if $\lambda$ represents the amount of saRNA (including activated and inactivated, where the proportion of activated saRNA in the total is $\theta$), $\mu$ represents the amount of protein, and $\varphi$ is a state function variable, then the joint distribution rate is:

the master equation of the system over time is:

where $\lambda, \ \mu \in \mathbb{N}_0$, and the influence terms in the differential equation are all influence terms added considering feedback regulation. When there is a large number of genes and proteins, according to the Law of Large Numbers, the result will converge according to the probability, but this requires a lot of differential equations to describe the phenomenon, as the number of molecules increases, the number of differential equations required will exceed the solvable range, so considering practical applications, it is necessary to simplify the above master equation and describe the system with stochastic differential equations to make the problem solvable[9].

For the above phenomenon, the state of the gene is described by the random process $G(t)$:

The stable solution of the above equation is described by expectation as well as variance, as shown in the figure:

To better describe the phenomenon, different molecular behaviors can be modeled separately over time:

3.2 Feasibility analysis based on deep learning

Artificial intelligence has been widely used in biological structure prediction, and the prediction of translated structure, especially the prediction of protein folding structure, has made great progress after the emergence of AlphaFold [10]. In this model, with the assistance of meta-learning, simple RNA sequence-based structure prediction is trained, and the feasibility of the project is analyzed through the results of structure prediction. Since AlphaFold does not solve the structural prediction caused by mutations, and mutations often occur in general biological experiments, the structural model based on mutation sequences is trained with the assistance of meta-learning, and the method of randomly adding restrictive white Gaussian noise is used in the generation of mutation raw samples, and the mutation sites are classified into different tasks in quantitative form, and the tasks are divided into Training Tasks and Testing Tasks in a ratio of 7:3. In each task, the original sample data is randomly produced by adding the method of restrictive white Gaussian noise for certain times, so that the data is more realistic, and the samples in each task are divided into Support Set and Query Set, according to which the Model Agnostic Meta Learning (MAML)[11] is used to train model parameters as the loss function to consider the value, and a random restart algorithm is added to avoid falling into local extremums in the optimization process. In the final protein structure prediction, this model uses the method of embedding AlphaFold to predict the results.

The criteria for evaluating the feasibility of the project are: after entering the genetic sequence in a certain plasmid measured in the experiment, the corresponding protein and its structure are predicted, and the experimental design is feasible if the final result (taking into account the bond angle, spacing, and amino acid type) matches the designed antibody with a high probability. By entering the sequence constructed in the experiment, the result feasibility is predicted.

The selection process of antibodies can be with the assistance of deep learning models, in the process of selecting the best antigen and antibody binding, the binding energy and chemical pairing relationship are synthesized, and finally the binding method with the highest score and the corresponding protein structure are selected (see the Deep Learning Section for details), this chapter only discusses the binding probability of antigens and antibodies under the premise of selecting antibodies. Antigen-antibody binding is essentially a ligand receptor binding problem, consider the following chemical equation:

The extent to which the reaction proceeds can be described by the dissociation constant:

Since the probability of antigen-antibody binding is occupied by the receptor rate, the probability of the complex occupying the whole is:

Considering the solution lattice model, i.e., $\Phi$ ligands are randomly distributed in $\Omega$ lattices, according to the Boltzmann distribution function:

The probabilities can be deduced as [6]:

where $\varepsilon_b$ is the ligand-receptor binding energy and $\varepsilon_{\text{sol}}$ is the energy of the ligand in solution, comparing (54)(56) two equations:

Taking $k_b = 1.38 \times 10^{-23} J/(K \cdot T)$, selecting the unit volume of $1nm^3$, $\varepsilon _b - \varepsilon_{\text{sol}} \approx -10 k_bT$ in the lattice model, it can be seen from the results that 95% binding rate can be achieved in a solution with a concentration of $10^1$ order of magnitude, and compared with different binding energy differences, the greater the binding energy difference, indicating that the antigen-antibody binding is more specific, the more stable the binding, and the lower the concentration required to achieve 95% binding rate. In this project, our antibodies have high specificity, and the reference concentration provided in the modeling can be selected to maximize antibody efficiency.