Advancing Plant Synthetic Biology Using Multiscale Mathematical Modeling

1. Introduction

The first green revolution successfully addressed the food crisis of the 1960s, by intensifying crop production under favorable environmental conditions with the aid of excessive fertilizers [1]. This approach led to an increase in the environmental impact of agriculture that will be charged to the current human generation [1]. These costs include land degradation with depleted soil’s nutrients, water eutrophication from fertilizer residues, over-extraction of groundwater, pest resistance caused by pesticide overuse, extinction of indigenous crop varieties, and increased greenhouse gases emissions that contribute to climate change [1,2,3]. Climate change triggers extreme weather conditions resulting in severe flooding and drought that jeopardize crop production in Asia [4]. Additionally, in sub-Saharan Africa, Striga hermonthica remains the most destructive pest, causing up to 100% yield loss in some crop fields [5]. In addition to the escalating climate crisis, ongoing conflicts, and recurring pest outbreaks are collectively undermining global crop production. For instance, the conflict between Russia and Ukraine is jeopardizing global maize production, as the two countries have a combined global production share of 17% for that cereal [6,7]. The continuous increase in world population also raises concerns about global food insecurity [8]. Projections by the Organization for Economic Co-operation and Development (OECD) and the Food and Agriculture Organization (FAO) of the United Nations (UN) estimate a world population of nearly 10 billion by 2050. This scenario is raising the stakes for sustainable food production and there are calls for a second Green Revolution that emphasizes production of crops that are resilient to suboptimal environmental conditions and create a more manageable environmental footprint [9].

To meet these challenges, researchers are turning to synthetic biology in plants, which employs biological engineering principles to improve beneficial plant traits in a way that goes beyond the speed and capacity of natural processes [10]. Through metabolic engineering and synthetic biology, plant scientists were able to remodel metabolic pathways to increase production of useful metabolites or inhibit synthesis of less useful ones [11]. Synthetic biology products that are now available in the market include Impossible Foods’ bleeding plant-based burgers, created from an engineered strain of the yeast Pichia pastoris that produces soy leghemoglobin and gives the burger a meaty sensation [12,13]. There is also Zymergen’s polyimide film, named hyaline, that is made from bio-sourced monomers [14], Pivot Bio’s biological nitrogen fertilizer for cereals [15], and Calyxt’s high-oleic oil from soybeans called Calyno [16,17]. More recently, printed vegan salmon also became available in supermarkets and many vegan fish alternatives are currently being developed [18]. As such, plant synthetic biology holds immense potential to steer agriculture towards a more sustainable future, alleviating global food insecurity.

While redesigning crops through plant synthetic biology is undeniably promising for the future of agriculture, it can also cause environmental damages. One of the risks that comes together with synthetic biology is producing a genetically engineered species that may increase ecosystemic problems once released into the natural environment [19,20]. Engineering plant genes may also initialize pleiotropic effects to the organism, due to off-target effects that can lead to unforeseen changes in the target species [19]. These risks can be mitigated using the combination of plant synthetic biology and mathematical modeling, which allows identification of mutational targets through the use of parameter estimation, bifurcation analysis, sensitivity analysis, and optimization [21,22].

To improve risk mitigation, integrated mathematical models should be applicable at different levels and with varying granularity. This is so because plants are complex systems composed of roughly modular components, and this modularity extends from the genome level to the tissue and organ levels [23]. However, the application of multiscale mathematical models, while frequent in microbial systems [24,25,26,27], remains underutilized in plants. Integrating technological advancement in the analysis of available agricultural data facilitates comprehensive understanding of crops, from genotype to phenotype to ecosystems, promoting effective environmental conservation and higher agricultural productivity [28,29].

To potentiate the creation and use of multilevel modeling in plants, this work reviews common modeling approaches and methodologies applicable to plants at various scales, from molecular to crop levels. We will focus on the most commonly used approaches, progressing from the simplest to the more mathematically complex methods. We also explore alternatives for integrating models across scales, discussing advantages and shortcomings of each approach. Our review concludes with insights into our ongoing efforts to create an integrated multiscale mathematical model for maize, modified through synthetic biology to be more resilient through enhanced resistance against Striga and droughts.

2. Model Building

Building multiscale mathematical models in plant synthetic biology remains a challenge, and as the model becomes more complex, so does its mathematical analysis. Further, integrating models of different levels and granularities poses additional complex challenges. In response to these challenges, we outline an iterative workflow for constructing a multiscale mathematical model in plant synthetic biology.

In Figure 1, we represent two nested cycles that summarize that workflow, portraying a minimalistic methodology for mathematical modeling that traverses across multiple scales from the molecular to the plant level. The flow in the inner cycle tackles the choice of methodology for model building and analysis, while the outer cycle portrays the integration between different scales of the plant processes one wants to model. The inner workflow applies individually at each step in the outer cycle to build a multiscale mathematical model of a plant. We should be able to use this model, once validated, to predict unintended pleiotropic effects of manipulating the genome on the plant at multiple levels.

In the remainder of this section we will review modelling strategies used in plants at various scales, from molecular to plant level. We will frame that review from the perspective of the workflow summarized in Figure 1. We note that the various stages proposed in Figure 1 may unfold to include more or less detailed modeling of the various plant subsystems. For example, shoot architecture may be modeled from a coarse level to a level where branching, leaves, and flowers are considered in detail.

2.1. Metabolic Modeling with FBA, Kinetic, or Hybrid Approaches

2.1.1. Applications

Metabolic modeling techniques, such as flux balance analysis (FBA), kinetic modeling, and hybrid approaches, play a crucial role in driving metabolic engineering towards efficiency [30]. Within the iterative design-build-test-learn (DBTL) cycle of metabolic engineering, metabolic modeling emerges as a valuable tool [11,31]. It offers a unique perspective on plant metabolism, yielding quantitative insights and informing new engineering strategies [31]. Particularly, it streamlines the trial-and-error cycles inherent in metabolic engineering, saving considerable time and resources [11,32,33].

We need to collect pathway data in order to build a metabolic model and mathematically represent it for analysis. When pathways are unavailable in databases, reconstruction and inference based on literature or genome analysis are necessary. This can involve utilizing gene sequences to identify related pathways and their data. Once we define the relevant pathways, we can associate a stoichiometric matrix to the system and use it for model construction. Additionally, preexisting kinetic information is often required for model creation and analysis. Table 1 provides a sample of databases containing kinetic and pathway data that are available to the reader for reference.

Most studies modeling plant metabolism employ Flux balance analysis (FBA)-based or kinetic-based models [31,33]. Flux balance analysis (FBA) stands out as a prevalent constraint-based approach for analyzing material and energy flow within metabolic networks [34]. Its advantage lies in the ability to analyze entire metabolic networks without detailed process mechanism information. FBA models have been extensively used to study plant metabolism, including in Arabidopsis thaliana, tomato, maize, and rice [35,36,37,38,39,40,41,42,43,44,45,46]. They facilitate understanding flux distributions and gene essentiality, requiring modest prior information.

However, FBA has limitations. It cannot accurately represent regulation and lacks information on transient responses or metabolite concentrations, crucial in metabolic engineering. For such cases, more detailed kinetic models are necessary, either developed from scratch or based on simplified FBA versions. Kinetic models offer mechanistic descriptions of metabolism, enabling regulation incorporation but demanding more prior information. They require knowledge of individual metabolic processes and parameter values for rate expressions, often obtained from databases like BRENDA, MetaCyc, or KEGG (Table 1). When using these databases, it is important to be aware that parameter values are frequently determined in vitro, and might be significantly different from in vivo values [33].

As FBA models, kinetic models can be used to analyze steady state flux distributions. In addition, they provide information about transient flux changes and concentrations, both at steady state and during transient responses. This is required in situations where metabolic engineers need to predict behavior of the biological system under genetic or environmental perturbations [33]. In addition, kinetic modelling enables sensitivity and stability analysis, which provide a framework for elucidating parameters responsible for the control of metabolic fluxes in a network [33,47,48,49].

Kinetic models typically focus on smaller parts of the metabolic network [31], employing rational expressions like Michaelis-Menten or Hill equations [50,51,52,53,54,55,56,57]. Other typical kinetic expressions include mass-action and Arrhenius equations [58], saturable and cooperative formalisms [49,59,60] as well as generalized mass action (GMA) kinetics [48,61]. GMA and other approximate representations are very useful when the mechanism for the individual processes is unknown [59,60]. Kinetic models were used to study, among other things, the flavonoid biosynthesis pathway in Arabidopsis thaliana [58] and tomato [54], central carbon metabolism in potato [55], central carbohydrate metabolism [53], aspartate-derived amino-acid pathways [51], and 2-C-methyl-D-erythritol 4-phosphate (MEP) pathway [56] in Arabidopsis thaliana, benzenoid network in Petunia hybrida [52], lignin biosynthesis in Brachypodium distachyon [61], and nitrate reductase activity [50].

FBA and kinetic models can sometimes be combined into hybrid models [33]. By integrating both methodologies, hybrid models overcome limitations and enhance model power. Creating a hybrid model is a two-step process that involves creating an FBA model to generate flux distributions, then using these distributions to estimate parameter values for kinetic sub-models. Hybrid models provide detailed concentration insights from kinetic modeling while validating results with FBA outputs, offering a more accurate representation of real-world metabolic processes [62]. Another parameter determination approach is ensemble modeling which calculates the parameter values with references to the steady state fluxes [31,63].

2.1.2. Creating a Metabolic Model

Figure 2 summarizes a typical workflow for metabolic modeling with FBA, kinetic, or hybrid approaches. Here, we illustrated the input, method/approach, and output of a plant metabolic model. Let us now look at the modeling process in more detail. In general, given that we a have a metabolic pathway with $m$ number of metabolites $M_i$ and n number of reactions/fluxes $v_i$, we can generate an $m$ by $n$ stoichiometric matrix $S$ [64,65,66,67]. After obtaining the stoichiometric matrix $S$, we can build a flux vector $\mathbf{v}$ with $n$ by 1 dimension where each row is the flux $v_i$. Then, we apply matrix multiplication to $S$ and $\mathbf{v}$ to get the metabolite/substrate vector $M$ with $m$ by 1 dimension where each row represents the rate of changes of each $M_i$. We represent this mathematically using Equation (1).

$${\displaystyle \frac{\mathrm{d}\mathrm{M}}{\mathrm{d}\mathrm{t}}}=\mathrm{S}\mathrm{v},$$

(1)

From equation (1), we diverge into different paths depending on the modeling approach. To create an FBA model, we assume a steady state of flux where the rate of changes in vector M is zero, that is ${\displaystyle \frac{dM}{dt}}=0$. Hence, Equation (1) becomes

$$\mathrm{S}\mathrm{v}=0.$$

(2)

At this stage we constraint the flux vector v using an objective function that could, for example depend upon the rate of material/energy entering and leaving the system and the minimal and maximal amounts of that flux that need to flow through different sub-processes for the system to be viable. We can then use linear optimization to solve the flux distribution and find the appropriate solutions to that optimization problem.

In constructing a kinetic model, we begin with Equation (1) and select suitable mathematical formalisms to represent the kinetics of individual processes and their dependence on substances directly influencing flux. Various formalisms, including linear, rational, and power laws, among others, are available [60]. Each formalism is associated with parameters whose values need to be estimated to predict the dynamic behavior of the model quantitatively. We can obtain initial parameter estimates for example from databases listed in Table 1 or from primary literature sources.

Equation (1) is equivalent to ${\displaystyle \frac{dM}{dt}}={\left[{\displaystyle \frac{d{M}_1}{dt}},{\displaystyle \frac{d{M}_2}{dt}},\dots ,{\displaystyle \frac{d{M}_m}{dt}}\right]}^T$ which is a nonlinear ODE. In analyzing nonlinear ODEs, a common practice is to approximate using linearization of the system. To linearize Equation (1), we calculate the Jacobian matrix (J) by getting the partial derivative of each ODE with respect to each dependent variable.

Using the Jacobian matrix (J), we determine the stability of the system through the computation of the eigenvalues of the matrix J. The real part of the eigenvalue of a stable system must be negative for the steady state to be stable. With this into consideration, we use different methods to understand how certain factors affect the concentrations of metabolites. For example, sensitivity analysis [65,66] helps us identify which parameters can cause significant changes, parameter scanning helps us gauge the extent of these changes, and time-series analysis helps us track how metabolite concentrations change over time [48,57].

From Equation (2), the output of FBA yields a flux distribution v that minimizes or maximizes the given objective function. Common applications of these models involve modifying constraints and initial conditions or altering the system by removing selected processes from the diagram, thus assessing their impact on the model's solution space. On the other hand, the output of the kinetic model encompasses all metabolic concentrations involved in the metabolic pathway, offering a comprehensive view of the system's behavior under specified conditions.

2.2. Shoot and Root System Architecture

Shoot and root plant system architectures are essential for resource acquisition during plant development. Yet, conducting plant experiments, especially in the field of plant synthetic biology, is inherently challenging with respect to time and cost. It is also hard to estimate the risk of disrupting local ecologies when introducing new plant varieties. Developing a 3D architecture visualization model is crucial to predict plant shoot/root growth and development, facilitating reasonable approximate estimations of crop productivity [68,69]. Thus, there has been a growing interest in building mathematical models of the shoot and root system architecture to do in silico investigation in plants. However, this field is still less mature than that of metabolic modeling and methodological development is ongoing.

2.2.1. Shoot System Architecture

2.2.1.1. Applications

The shoot system architecture (SSA) is generally the above ground portion of the plant that consists of stem, leaves, buds, flowers, and fruits. Optimizing the shoot growth of wheat revolutionized agriculture during the 1960s resulting in the First Green Revolution. First, the agronomist Norman Borlaug developed a high yield and disease resistant wheat variety [70,71]. However, the wheat stalk of this variety was too thin to support the weight of its own grain, leading to a type of breaking called lodging. Thus, he further developed a dwarf variety of this wheat which can bear heavier weight [71]. This development of dwarf wheat solved an impending global food crisis during that time.

Creating dwarf wheat through traditional breeding took many years [72]. With the current gene editing technology, it is already possible to speed up the process of developing more desirable plant varieties, creating specific modifications in the genome of plants to obtain specific traits. Using models to simulate the effects of these genetic modifications on the phenotype of SSA could mitigate risk and predict the quantitative effect of genetic interventions on shoot traits and crop yield more accurately. This could speed up the process of producing new crop varieties.

There are several models for SSA that have been developed over the years [73,74,75,76]. L-systems underlie many of these models [73,76,77,78,79,80,81,82,83,84,85,86,87]. This modeling formalism was named after Aristid Lindenmayer, and it is widely used to model plant morphogenesis [88]. In Table 2, we listed several functional-structural plant modeling (FSPM) platforms from the Quantitative plant website [89]. We can also use these platforms to model root and whole plant architectures.

2.2.1.2. Creating a Shoot System Architecture

In Figure 3, we present a workflow for building a SSA model. Initially, we need to determine whether the plant species is monocot or dicot, as they have different architecture from one another. A monocot generally has tillers, with secondary shoots branching from its primary shoot. A dicot typically has only one primary shoot [90].

In modeling SSA, we need to describe the phytomers, which are the fundamental structural units of the shoot. These phytomers are produced from shoot apical meristems (SAM), also called terminal bud [74,76,91]. Figure 3 shows an example of phytomers. There, each phytomer is a stem segment that is composed of a leaf, an axillary bud, a terminal bud, a node, and an internode. The phytomers are stacked up to form a shoot, whose form may vary depending on the plant species. Wen et al. [76] defined a geometrical model of a shoot as

$$Shoo{t}_i^n=\bigcup _{j=1}^n{R}_{ij}·T_{ij}·Phytome{r}_{ij}.$$

(3)

Here, $R_{ij}$ and $T_{ij}$ describes the rotation and translation of the ${Phytomer}_{ij}$ and the $Shoo{t}_i^n$ is the union of $n$ number of phytomers.

We can use these models to quantify shoot architecture traits such as leaf area (an important factor for capturing light), plant height, stem length, volume, etc.

2.2.2. Root System Architecture

2.2.2.1. Applications

The most frequently used root phenotyping techniques include shovelomics, the monolith method, and X-ray computer tomography [92,93,94]. Shovelomics [92] is a rapid and high throughput phenotyping method. This is done by digging a shovel size portion of the root and then manually measuring 10 root architectural traits. The monolith method [93] uses a steel cylinder with 20 cm diameter and 30 cm depth for excavating roots and measuring their phenotypes. This method allows excavation of deeper roots and uses a root image analysis software to quantify root traits. More recently, X-ray computer tomography (CT) was used for 3D visualization that captures more detailed root architecture traits [94].

In comparison to shoot, root data gathering is particularly challenging. This difficulty exponentially increases as the roots continuously grow deeper into the ground. Thus, the development of easy-to-use methodologies to facilitate root phenotyping is still an important unmet need. The belowground disposition of the root makes it more beneficial to simulate a 3D visualization of the root system architecture (RSA), when compared to SSA. Similarly to SSA, the most common approach in modeling RSA is also the L-system [88,95].

Researchers have long been engaged in constructing virtual models of roots and on using those models to simulate the effect that changing environmental conditions has on the RSA, with the aim of prioritizing interventions that might increase crop yield [96]. The development of root modeling platforms dates back to the 1970s [97,98]. Later, more advanced modelling platforms were created. ROOTMAP [99] is a notable 3D model focusing on fibrous root systems and their proliferation. Other platforms, like RootTyp [100,101] and DigR [102], highlight the influence of soil in root development, while SPACSYS [103] delves into root-soil interactions and their impact on crop yield.

Unlike metabolic models or SSA models, many RSA modeling platforms are not publicly available [104]. This hinders the advancement of the field. The rise of open science spurred initiatives like OpenSimRoot [105] , an open-source 3D RSA model accommodating various plant species and focusing on nutrient acquisition's feedback to root development. Similarly, CRootBox [106], derived from RootBox [95], fills the gap between in situ and in silico studies by integrating 3D RSA simulations into experimental setups such as cylindrical pots and rhizotrons.

These platforms contribute to untangling the complexities of RSA. More recently, a modeling study [107] proposed an approach for predicting root/shoot system architectures. Their iterative and stochastic method generates elongation length, root direction, and branching until building an individual plant root architecture. Notably, their approach allows for modular introduction of hormonal and metabolic factors influencing architecture development, predicting their effects on various root traits like length, branching, volume, and surface area.

2.2.2.2. Creating a Root System Architecture

Figure 3 summarizes a three-step workflow for building models for root system architectures (RSA). In a parallel situation to that of plant shoots, plant roots may be classified either as tap root for dicots or fibrous roots for monocots. Taproots have one primary root (or main root axis) while fibrous roots have several additional root axes, for example seminal and nodal roots. As such, the initial step in creating a model for a RSA is to decide which type of root one needs to model.

Subsequently, we need to define a root axis for the RSA model. Figure 3 illustrates a root segment of the RSA, consisting of a main root axis that displays secondary root branching in the form of lateral roots. This main root axis accurately represents a taproot. For modeling fibrous roots, we need to generate additional axes that describe seminal and nodal roots.

A step-by-step description of such a model building process can be found in a recent literature [107]. That process requires defining at least two parameters: the growth elongation rate and the growth angle, which are important for growth elongation rate and growth angle, respectively. We initialize the root system architecture is initialized at the origin in 3D space. Iteratively, the approach decides the magnitude of the elongation and direction of the root for each step. Stochastically, we define branching points from where secondary roots will branch and grow, generating a full architecture. This model outputs a 3D visualization of the RSA of a plant, something hard to image experimentally. After generating the RSA, we can also quantify length, branching, volume, surface area, and other traits of the RSA.

2.3. Resource Acquisition

Resource acquisition in plants encompasses obtaining energy and fixing carbon via photosynthesis and collecting mineral and fixing nitrogen through the roots. Modelling plant resource acquisition is often but not always connected to modelling RSA and SSA. Understanding how changes in SSA or RSA influence light acquisition via photosynthesis and water and mineral nutrient uptake from the soil is usually the primary objective of resource acquisition models.

In the context of resource acquisition, researchers use 3D SSA/RSA models to identify factors that change the architectures in ways that maximize a plant’s resource acquisition. Examples include modelling the effect of changing leaf’s area on the amount of light captured by the plant [108,109,110], or simulating the effect that changing spread and depth of root into the soil has on mineral, nitrogen, or water uptake through the root [111].

2.3.1. Photosynthesis

2.3.1.1. Applications

Photosynthesis sustains life on our planet through the supply of food and oxygen [112], allowing plants to grow, develop and reproduce [113]. This process lets plants capture energy from light, use that energy to fix atmospheric carbon from CO₂ into sugar, and release oxygen into the air as a byproduct.

Photosynthesis has two phases. The first phase converts energy from light into a chemically usable form. The second phase fixes atmospheric CO₂. By and large there are three different types of metabolism to fix that CO₂: $C_3$-, $C_4$-, or Crassulacean acid metabolism (CAM) [113,114,115].

Figure 4 illustrates how to model photosynthesis. The process initiates with the decision regarding which type of metabolism needs to be modeled ($C_3$–, $C_4$ –, or CAM). Model inputs always include (but are not limited to) light and $C{O}_2$. Regardless of the type of photosynthetic metabolism, three broad approaches can be used to build the model: empirical, mechanistic, or hybrid approaches that combine the first two [113]. In addition, and independent of the approach, the models can be either deterministic or stochastic and simulations can run on continuous or discrete time steps. This study [113] highlighted that decisions regarding which approaches to use should consider minimizing the number of model parameters and correctly representing physical, chemical, and biological laws. One should also consider minimizing the variance of the simulation outputs and the deviation between predicted model values and measured experimental values.

2.3.1.2. Creating a Photosynthesis Model

The Farquhar – von Caemmerer – Berry (FvCB) model pioneers mathematical modeling of photosynthesis of $C_3$ plants [114]. The FvCB model calculates the $C{O}_2$ assimilation using Eq 4 [113,114,116]:

$$A=V_C-0.5V_O-R_d.$$

(4)

Here, $V_C$ and $V_O$ represent carboxylation rate and oxygenation rate, respectively, and $R_d$ is the mitochondrial respiration rate.

C₄ photosynthesis relies on the coordination between mesophyll and bundle sheath cells within leaves. Initially, CO₂ is fixed into a C₄ acid, likely in the mesophyll, by phosphoenolpyruvate carboxylase (PEPC) [115,117]. This 4-carbon compound then moves to the bundle sheath, where it undergoes decarboxylation, providing CO₂ for Rubisco to initiate the regular C₃ carbon fixation process. Therefore, modeling C₄ photosynthesis entails additional initial steps, along with a mathematical model of C₃ photosynthesis (see [115] for details).

Plants in arid regions, characterized by less than 25 cm of rainfall per year, have adapted to dry environments by evolving CAM metabolism, a type of photosynthesis that fixes CO2 more efficiently while reducing evapotranspiration. CAM photosynthesis operates through a circadian metabolic cycle, with stomata opening at night and closing during the day to minimize water loss. During the night, PEPC fixes CO2 to produce malic acid, which is stored in the cell vacuole. During the day, the vacuole releases malic acid, which is then decarboxylated through the Calvin Cycle. Recent review [118] presented several mathematical models of CAM photosynthesis, including ODE-based models [119,120,121] and FBA – based models [122].

In photosynthesis models, the typical output variable is the rate at which plants fix atmospheric CO₂ and/or the rate at which plants use that CO₂ to produce sugars and other compounds. Various definitions of photosynthetic rate exist in the literature, all linked to the plant's use of atmospheric CO₂ to produce sugars [113,114].

2.3.2. Nutrient Uptake Model

2.3.2.1. Applications

One of the primary functions of roots is to acquire water and nutrients from the ground and supply them to meet the metabolic needs of plants. Nutrient acquisition models are normally layered on top of RSA in silico models [105,106,123,124]. These models are very useful to investigate the effect of rapidly changing environmental conditions on plant growth, and prioritize which costly field trials to make [125], saving time and money. They are also useful in combination with traditional approaches to develop innovative practices in fertilizer management [126]. In addition, quantifying root nutrient uptake provides insights into potential crop yield.

Over the years, the development and application of several models to quantify root nutrient uptake has been on the rise, aiming to optimize fertilizer usage in agriculture. Initially, nutrient uptake was modeled as directly proportional to the available nutrient concentration in the soil, with root absorption power serving as the constant parameter [127]. However, solving the system of partial differential equations (PDEs), whether analytically or numerically, posed a significant challenge. To address this issue, a study [128] introduced the Crank-Nicolson method to numerically solve their nutrient uptake model. Building upon this approach, another study [129,130] incorporated root system growth into their models using the same numerical method. Later, the roots of rapeseed were modeled to investigate the influence of root hairs on the uptake of phosphorus [131]. These studies were extended [132] to other plant species, who quantified how plant phosphorus uptake strongly influenced the number, length, and radius of root hairs.

The Nye-Tinker-Barber (NTB) model became the most frequently used for modeling plant nutrient uptake [133,134]. A study [135] extended this model, demonstrating the existence of an explicit closed-form solution to the diffusion-adsorption problem [136,137]. Additionally, it was further extended [138] and analytical solutions were provided to convection-diffusion equations applicable to general solute nutrient uptake scenarios.

Recent advancements include the use of perturbation expansion methods to approximate nutrient flux and concentration solutions [139]. This approach divided the rhizosphere into inner and outer fields, aligning with the root surface, and derived approximate analytical solutions for nutrient uptake flux at the root surface and global nutrient diffusion of the NTB model. Another study [140] contributed to finding analytical solutions using Laplace transforms.

2.3.2.2. Creating a Root Nutrient Uptake Model

Figure 4 illustrates the modeling process of nutrient uptake, beginning with the definition of the biological problem. This involves identifying the specific nutrient being modeled, which can be categorized as either mobile (such as nitrogen, sulfur, boron, manganese, and chlorine) or immobile (including phosphorus, potassium, magnesium, calcium, copper, iron, and zinc) in soil. A classification of the mobility of mineral nutrients in soil [141] is presented in Table 3.

Following the classification of nutrients, diffusion equations mathematically represent the transport process from the soil to the root. These equations are typically described using PDEs, which incorporate both time and space as independent variables [142]:

$${\displaystyle \frac{\partial c\left(x,t\right)}{\partial t}}={\displaystyle \frac{\partial \mathit{J}\left(x,t\right)}{\partial x}},$$

(5)

where the concentration (c) changes over time (t) and the flux (J) depends on the spatial variable (x).

This simple diffusion equation becomes more complex as additional spatial features of solute transport are incorporated. A review study [143] categorized various solute uptake models from the literature based on their complexity, including those representing diffusion only, diffusion with advection, diffusion with reaction, and diffusion with both advection and reaction.

Once we have mathematically represented root nutrient uptake using PDEs, the next step is to determine their solution, which can be either analytical or numerical. While finding analytical solutions can be challenging, a widely used method is non-dimensionalization, which simplifies PDEs by introducing non-dimensional variables [135]. Non-dimensionalization facilitates the process of finding an analytic solution [144]. When analytical solutions are hard to find, we use simulations to find numerical solutions, often by discretizing time and space. A study [145] extensively reviewed different numerical methods for discretization.

After determining the solution, various analyses can be performed, including parameter estimation, sensitivity analysis, or prediction of soil nutrient uptake fluxes and concentrations.

2.4. Shoot-Root Interaction

2.4.1. Applications

The allocation of resources in plants is a crucial ability that determines their survival in challenging environmental conditions. Current climate trends emphasize the need for crops that have phenotypic plasticity and can adapt to growing with less water and more heat. This phenotypic plasticity comes from the capacity of a given genotype to express different phenotypes in response to varying environmental conditions [146]. Modeling shoot-root interactions is especially useful for simulating different scenarios and predicting the range of phenotypic plasticity of plants in different surroundings. In addition, modeling shoot-root dynamics allows for further exploration of plant productivity, considering the allocation of nutrient minerals from uptake and carbon from photosynthesis. This integration creates a multiscale mathematical model capable of facilitating in silico experiments ranging from the metabolic level to the whole plant level.

Researchers proposed many transport-resistance (TR) models to analyze plant growth dynamics. A TR model that pioneered this field [147] introduced a leaf-stem-root model that accounted for the utilization and transportation of photosynthate among different parts of the plant. In a subsequent study [148], they constructed a shoot-root model focusing on the exchange of photosynthetic output between the shoot (carbon) and root (nitrogen) uptake. Additionally, he incorporated the chemical and biochemical conversion of these substrates into structural dry matter [149]. This was expanded [150] by considering the effect of water potential and transpiration in his model. He utilized the Münch flow mechanism of phloem translocation and introduced a switch based on concentration gradients. In a more recent study focused on phosphorus dynamics [151], researchers combined experiments and models to investigate shoot-root interaction over a day-and-night cycle, using Petunia hybrida as a model organism.

There are other whole plant level models that considers root-shoot interaction and are not TR models. For example, CPlantBox [152] and PiafMunch [124,153] applied cohesion-tension theory for water flow and Münch flow theory for carbon flow.

2.4.2. Creating a Shoot-Interaction Model

Figure 5 outlines the modeling process for shoot-root interactions, depicting the input compartment consisting of resources acquired from photosynthesis and nutrient minerals from root uptake. These models primarily focus on resource allocation, emphasizing the exchange and transfer of resources between the shoot and root. In describing the utilization and transportation of substrate, most models employ the TR approach. This modeling approach is favored for its realistic representation of transport, accounting for substrate losses during transport and translating C/N fixation into root growth, using functions that convert volume into plant structure [154].

Still, different models use alternative strategies. For example, the teleonomic model [148] prioritizes factors like the shoot-root ratio, growth, and substrate concentrations based on environmental conditions [155]. Other models give more importance to other factors, such as soil water potential and root temperature effects on plant partitioning of photosynthate [156].

Regardless of the approach used, shoot-root models need to account for the interactions between two interconnected compartments: the shoot system and the root system [151]. The shoot part may range from simple shoot concentrations of substrates to more complex models specifying shoot concentrations of stored starch and soluble sugars, including shoot structure/architecture and mineral nutrients from the roots such as nitrogen, phosphorus, and potassium [124,152,157]. Similarly, the root part comprises exchanged substrate concentrations from the shoot, root structure/architecture, and nutrient uptake from the soil [151,152]. The exchange of soluble sugars from the shoot and mineral nutrients from the root is driven by a transport function [149], with the growth rate dependent on the substrate concentrations it contains.

These models often output the concentration of different substrates for shoot and root over time, mathematically represented by ordinary differential equations (ODEs). Thus, they enable in silico stability and sensitivity analysis of the shoot-root ratio under various environmental conditions.

3. Discussion

3.1. The Potential of Multiscale Mathematical Modeling

Plant synthetic biology aims to engineer plant biological systems to improve crops through direct manipulation of genetic circuits in plant genomes [10]. Given the complexities of such manipulations, they come associated with potentially serious biosafety risks [19,158]. Multiscale mathematical models can be a valuable tool for risk mitigation through the in silico simulation and analysis of those effects. These models provide another layer of quality control for the design stage of design-build-test-learning (DBTL) cycle of synthetic biology [159]. Ultimately, this approach can expedite the development of new crop varieties and reduce the costs associated to that development [160].

While crucial to test the effect of interventions, traditional plant breeding techniques and experimental studies are time-consuming and unpredictable [71]. We can use multiscale mathematical models to predict the effects of genotype modifications on the plant’s phenotypic traits. Integrating in silico investigation through multiscale mathematical models into the DBTL cycle of synthetic biology enhances speed and accuracy of predictions made about for the effect of genome manipulation on in planta implementation [160,161]. In addition, we can fine-tune their precision in quantifying the changes caused by the genetic intervention to molecular, architecture, resources, and whole-plant levels by incorporating experimental information that becomes available during the model building process.

Understanding and correctly representing the interplay between those levels is crucial for improving plant growth and productivity and optimizing fertilizer management while reducing environmental impacts. For instance, growing shoots and roots are responsible for resource acquisition, which is then used to produce sugar that is consumed for the shoot and the root to continue to grow. Several studies [162,163,164] pinpoint the influence of carbon, nitrate, and sugar in shoot and root development. This and many other nonlinear internal feedback loops naturally exist in biological systems [57,165]. Accurately predicting the effect of those regulatory loops on the dynamic behavior of the plant requires using multiscale mathematical models. The predictions can subsequently be validated in comparison to experimental evidence [28]. There is a review [166] of several designs of synthetic gene circuits for plants which can be modeled using a multiscale mathematical model [28] and following the methodology that we presented in Figure 1.

Account for the complexity of plant systems requires multiscale mathematical models, as single-scale models cannot capture the intricate dynamics of the whole system [167]. Figure 1 illustrates an iterative process for creating these models, allowing entry at any step of the model building process. This approach helps understand how molecular changes propagate throughout the organism, informing strategies for genome modification to achieve specific objectives.

Despite the beneficial contribution of this approach to the advancement of plant synthetic biology, there is a significant gap between scales when it comes to available knowledge and methodological tools for creating and analyzing models at each scale. This is evident when comparing for example metabolic level models to shoot/root system models. Such gaps complicate the integration between scales within a model. This integration is crucial for evaluating how genetic manipulations interact with environmental and metabolic cues, influencing plant development. The modular nature of plant systems suggests an approach to facilitate scale integration in multilevel models: make the models representing the different levels as modules [27,28,168]. Also, there is a continuous effort to build Crops In Silico [167], which is a platform for plant mathematical models highlighting the value of multiscale approach.

A critical aspect in doing so is describing the connections between scales or translating the output of a model from one scale as input into another. One solution is to represent each scale in a separate submodel, using appropriate mathematical representations and tools, and interconnect them through relevant variables. This approach has been successfully applied to create multiscale mathematical models of several microorganisms [24,25,26,27].

The development of useful multiscale plant models depends on accurate parameter measurements. Advanced plant phenotyping technologies [92,93,94] are leading to an era of big data, providing modelers with highly accurate measurements. A study [169] proposed a modeling workflow centered on digital twins, using X-ray computed tomography (CT) to capture detailed root images for constructing root skeletons and RSA models. This approach allows for comprehensive analysis of plant phenotyping data, providing deeper insights into various plant functions [170]. While another study [107], modeled of the effects of varying concentrations of strigolactones on RSA demonstrated changes in root elongation and lateral root branching density.

The construction of whole-plant digital twins is likely to be forthcoming. Building whole-plant digital twins will likely be facilitated by artificial intelligence (AI). AI algorithms can analyze experimental data to inform model-building choices and later analyze data generated from digital twins, offering insights into plant growth, development, and responses to environmental stimuli. Leveraging digital twin technology alongside AI has the potential to optimize agricultural processes, improve crop yields, and contribute to sustainable food production. This holistic approach enhances our ability to analyze complex plant systems and paves the way for transformative developments in agriculture, addressing global food security challenges.

3.2. Limitations of Multiscale Mathematical Modeling

Modeling complex plant systems requires numerous parameter values, which can constrain the predictive power of models. Technological limitations often make it challenging to obtain certain parameter values, particularly at the molecular level. Detailed multiscale mathematical models also require substantial computational resources, leading to time constraints. As technology evolves, overcoming these limitations will be pivotal for realizing the full potential of multiscale mathematical modeling in plant research.

The integration of different scales within a multiscale framework is essential but challenging. Generalization of submodels is possible, but requires adaptation of the models in ways that depend on the study's specific objectives. For instance, the whole-cell model [24] offers intricate detail compared to a strigolactone biosynthetic model [57], reflecting distinct research goals.

A significant barrier in plant modeling is the disconnection between wet lab and dry lab researchers. Experimental findings often do not provide the necessary parameter values for mathematical models. Bridging this gap requires collaboration between researchers from both domains, facilitating data-driven model building and model-assisted experiments. This interdisciplinary approach is crucial for advancing plant modeling.

3.3. Concluding Remarks

Multiscale mathematical modeling holds immense potential for characterizing biological phenomena and enabling targeted bioengineering for improved crop productivity and sustainability [10,28,167]. However, integrating single-scale models into a comprehensive multiscale framework remains a critical challenge. Exploring diverse modeling approaches at each scale and harmonizing them into a unified whole-plant model is essential.

Deciphering and integrating different single-scale models are ongoing challenges in constructing a multiscale mathematical model. This review proposes a framework for constructing individual submodels tailored to plants within plant synthetic biology. The application of multiscale mathematical modeling in plant biology is a fertile area for future exploration. It has already proven invaluable in enhancing our understanding of plant complexity and driving advancements in agricultural research. Encouraging more applications of multiscale mathematical modeling in plants can contribute to productive and sustainable agriculture.

In summary, integrating multiscale mathematical modeling approaches promises to revolutionize our understanding of plant biology and enhance agricultural productivity. By leveraging diverse modeling techniques and frameworks, researchers can gain deeper insights into the intricate dynamics of plant systems, leading to more informed decision-making and innovative solutions for global food security challenges. Realizing this potential will require interdisciplinary collaboration, technological advancements, and ongoing refinement of modeling methodologies to effectively address the complexities of plant biology.