Tutorial 1: Package basics
In this tutorial, we'll cover the basics of this package. The main topics are:
Introduction: the concept of the package.
Modeling context
Creating a modular plant model
Defining structural and functional plant modules
Specifying model parameters
Coupling structure and function
Running the model
Visualizing the results
Setup
Before we can get started, we need to do some basic setup: activating our project and loading in the required packages.
using Pkg; Pkg.activate("../..")
using PlantModules
using PlantGraphs, ModelingToolkit, OrdinaryDiffEq, Plots
using PlutoUI; TableOfContents()
Introduction: modeling with PlantModules.jl
As the name suggests, one of the most important goals of PlantModules is to enable the user to easily model plant growth in a modular manner. In order to achieve this, we will define our plant in function of a few sets of similarly behaving parts or "modules" for short. Parts of the same module can be similar on the structural - and on the functional level.
On the structural level, some examples are:
An oak tree can be considered a repeating module in a forest.
A branch can be considered a repeating module in a tree.
A collenchyma cell can be considered a repeating module in a branch.
On the functional level, on the other hand:
Most structural modules of a tree (branches, leaves, root segments, etc.) contain water which will flow between them as dictated by hydraulic laws, which plays an important role in plant growth. All these structural modules share this same functional module.
One or more of the plant's structural modules will assimilate carbon through photosynthesis, while the others do not. On the cellular level, only cells containing chloroplasts should share this functional module.
As such, the workflow to create a model in PlantModules boils down to defining these modules and how they interact.
Example problem definition
For the first tutorial, we will model the response of a hypothetical tomato plant to different transpiration rates. In higher plants, an increase in transpiration rate causes a reduction in the water potential of the stem. This phenomenon can be explained by the cohesion-tension theory, which states that the water in plant vessels forms a continuous network through its strong cohesive forces, and pressure reductions caused by leaf transpiration cause pressure gradients throughout the entire network. We will verify here that our basic functional modules can capture this phenomenon.
Defining the model
Structure
The plant
The first step in modeling our plant's growth is defining its structure. For the first tutorial, we'll consider a very simple example: a tine tomato plant with a small stem and two leaves.
PlantModules uses graphs to define the relationships between modules. As such, structural modules need to be implemented as graph nodes of some sort. We'll be using the graph implementation from the PlantGraphs package for this tutorial, though any graph implementation can be used.
For our example, we will define two structural modules: a Stem module and a Leaf module. We will give both modules an attribute D, which corresponds to the name of the variable expressing the dimensions of the structural module. This is a vector containing one to three numbers, depending on whether the plant part is modelled as a sphere, a cylinder or a cuboid. The Leaves will be represented as flat cuboids, while the stem will be represented as a cylinder.
struct Stem <: Node end
struct Leaf <: Node
D::Vector
end
We now define the plant structure using a graph. By specifying the value of the attribute D, we can specify different initial sizes for the different plant parts. Note that the cylindrical stem only needs 2 values specified (radius and length), while the cuboid leaves need 3 values (length, width and thickness).
plant_graph = Stem() + Stem() +
(Leaf([5.0, 3.0, 0.05]), Leaf([1.0, 0.5, 0.05]));
plotstructure(plant_graph)
The environment
Plants need an environment to grow in. For most plants, the most basic environmental compartments that need defining are the soil, from which the plant gets water and nutrients, and the air, with which the plant exchanges gasses. We can make these compartments as complex as we want, such as different soil compartments for different layers, but we will stick to one graph node for each for now.
struct Soil <: Node end
struct Air <: Node end
soil_graph = Soil();
air_graph = Air();
We also need to define the connections between these environmental modules and the plant. Sadly, connecting both leaf nodes to the air would introduce a cycle into our graph, which is something many plant-focused graph implementations do not support.
To overcome this, PlantModules allows the user to input multiple graphs and then specify how to connect them:
graphs = [plant_graph, soil_graph, air_graph];
intergraph_connections = [[1, 2] => (getnodes(plant_graph)[1], :Soil), [1, 3] => (:Leaf, :Air)];
The above syntax means: "connect the first node from graphs[1] to all Soil nodes from graphs[2]", and all Leaf nodes from graphs[1] to all Air nodes from graphs[3]. Note that the order matters! For example, [1, 2] => (:Root, :Soil) is not the same as [2, 1] => (:Root, :Soil).
The full structural connection information can then be acquired by bundling the graphs and connections into the PlantStructure wrapper.
plantstructure = PlantStructure(graphs, intergraph_connections);
plotstructure(plantstructure)
For our tiny example plant, we only have one structural module that actually repeats, somewhat defeating the purpose of modelling in a modular manner: we may as well write this entire model out by hand! Rest assured, however, that the approach we're seeing here also works for larger plants, as we will see in the next tutorial.
Function
Now that we have an idea of the plant's structural modules, we need to assign them some kind of functionality. PlantModules expresses these as sets of differential equations, implemented in ModelingToolkit.jl. Assigning functionality then comes down to writing out the desired functional processes as ModelingToolkit models (or choosing them out from the prebuilt ones), mapping them to the structural modules, and finally specifying the initial- and parameter values for the different components of the system.
Defining functional modules
There are two types of functional module: next to the typical functional modules describing the function of a structural module, there are also the connection modules that describe how different structural modules are connected. Both are defined as ModelingToolkit.jl models. For this tutorial, we will limit ourselves to using functional modules already implemented in PlantModules.jl; implementation of new processes will be discussed in the second tutorial.
We list the pre-implemented functional modules here with a short description of each.
Hydraulics
hydraulic_module: describes the hydraulics-driven growth of a compartment.environmental_module: describes a compartment that is subject to hydraulic laws but does therefor not change in size.
Carbon dynamics
constant_carbon_module: describes a compartment with a constant amount of carbon.simple_photosynthesis_module: describes a compartment with carbon being produced during day and consumed proportionally to its concentration.
Total water potential
Ψ_soil_module: describes an empirical relationship between the total water potential and relative water content of soil.Ψ_air_module: describes the relationship between the total water potential and water potential of air.
Hydraulic conductivity
K_module: describes hydraulic conductivity as proportional to the cross area of the compartment.constant_K_module: describes a constant hydraulic conductivity.
Connection modules
hydraulic_connection: describes a water flow between two compartments based on the specific hydraulic conductivities of those compartments.daynight_hydraulic_connection: describes a water flow between two compartments based on the hydraulic conductivities of those compartments that decreases during night.constant_hydraulic_connection: describes a water flow between two compartments and specifies its own hydraulic conductivity.
Please note that the carbon dynamics modules are very simplistic and only exist so users can make a model that works out of the box. Users that want to create a realistic model are expected to implement some functional modules of their own in accordance with their own use for the package, as will be discussed next tutorial.
Mapping function to structure
When the functional modules have been defined, we can assign them to the different types of structural module and the connections between them. We can do this for both with a simple Dictionary. Note that, for now, while structural modules can be assigned multiple functional modules, their connections can only be assigned one connection module.
module_coupling = Dict(
:Stem => [hydraulic_module, constant_carbon_module, K_module],
:Leaf => [hydraulic_module, simple_photosynthesis_module, K_module],
:Soil => [environmental_module, Ψ_soil_module, constant_K_module],
:Air => [environmental_module, Ψ_air_module, constant_K_module]
);
connecting_modules = Dict(
(:Soil, :Stem) => hydraulic_connection,
(:Stem, :Stem) => hydraulic_connection,
(:Stem, :Leaf) => hydraulic_connection,
(:Leaf, :Air) => daynight_hydraulic_connection,
);
plantcoupling = PlantCoupling(; module_coupling, connecting_modules);
Defining parameter and initial values
Next to the equations themselves, an important part in describing the plant's functional processes is defining the equations' parameter - and initial values. Since we'll restrict ourselves to hydraulic equations for this tutorial, these parameters and variables are the ones described in the section Theoretical overview.
In PlantModules, parameter - and initial values are defined in a hierarchical manner, here listed in ascending priority:
Model-wide default values < module-specific default values < node-specific values
For our current example, this means that there is no point in specifying the initial values for the compartment dimensions D for our Leaf compartments in the module-specific default values or the model-wide default values, since we already defined these values in all of the Leaf nodes of the graph, which have the highest priority.
Before we change any of the default parameter - and initial values, we can take a look at the defaults used by PlantModules:
PlantModules.default_values
Dict{Symbol, Any} with 18 entries:
:η_night => 0.1
:D => [0.5, 5.0]
:h => 0.0
:M_c => 0.05
:W_r => 0.8
:K => 1000.0
:A_max => 2.0e-6
:W_max => 1.0e6
:M => 0.0003
:K_s => 10.0
:E_D => 50.0
:Γ => 0.3
:t_sunrise => 8
⋮ => ⋮
Now imagine we have some information on our plant in question and it calls for different default values. These changes will be inputted as a Dictionary:
default_changes = Dict(:Ψ => -0.2);
Next to changing functional values over the entire model, some of the information we have may require different default values specifically for certain structural modules. We can specify this as shown below.
module_defaults = Dict(
:Stem => Dict(:D => [0.5, 10.0], :M => 200e-6),
:Leaf => Dict(:shape => Cuboid(), :M => 300e-6),
:Soil => Dict(:W_max => 1e3, :K => 1.0),
:Air => Dict(:W_r => 0.7, :K => 1e-3)
);
As you can see, we changed the default \(shape\) for the leaf modules, as well as the initial dimensions \(D\) and metabolite concentration \(M\) of the stem. For both environmental compartments, we changed the hydraulic conductivity, and additionally the maximum water content for the soil comparment. For a detailed explanation of the parameters and initial values here, we again refer to the theoretical overview.
Values for the connection modules can be modified in a similar way:
connection_values = Dict(
(:Leaf, :Air) => Dict(:t_sunrise => 6, :t_sunset => 22),
);
In contrast to the nodes, PlantModules currently does not support edges of the same type (that is, between the same type of structural module) to specify different parameter - or initial values in the graph definition.
Finally, all parameter changes are wrapped in the PlantParameters struct.
plantparams = PlantParameters(; default_changes, module_defaults, connection_values);
Bringing it all together
Now that we have all parts of our model defined, all that's left to do is putting it together. For this we use the main workhorse of PlantModules: generate_system.
system = generate_system(plantstructure, plantcoupling, plantparams);
This function will generate the ODESystem describing the model. It is possible to fine-tune the model even further at this stage as described in the Customizing the model section of the docs, thought this should generally not be required.
Running the model
Running the model is taken care of DifferentialEquations.jl. Users that are unfamiliar with the package may want to take a brief look at the DifferentialEquations.jl docs.
time_span = (0.0, 48.0);
prob = ODEProblem(system, [], time_span, sparse = true);
sol = solve(prob);
In order to compare the results of higher transpiration, we remake and solve the problem with lower relative air humidity. This can be done by repeating the parameter, problem, system and solution definition with a different value for W_r, or we can simply use ModelingToolkit.jl's remake function. To make it easier to get the correct variable from our system with many subsystems, the get_subsystem_variables convenience function exists that extract all variables corresponding to a certain name and structural module (or connection between two).
air_W_r = get_subsystem_variables(system, plantstructure, :W_r, :Air)[1]
$$\begin{equation} \mathtt{Air6.W\_r}\left( t \right) \end{equation}$$
prob2 = remake(prob, u0 = [air_W_r => 0.5]);
sol2 = solve(prob2);
Showing the results
First we verify that there is indeed a higher transpiration in our second simulation by plotting the total incoming water influx of the air.
plot(
plotgraph(sol, plantstructure, varname = :ΣF, structmod = :Air, title = "Low transpiration", xaxis = "t (hr)"),
plotgraph(sol2, plantstructure, varname = :ΣF, structmod = :Air, title = "High transpiration", xaxis = "t (hr)"), ylims = (0.0, 0.25), yaxis = "ΣF"
)
Finally, we can inspect whether the effect on the water potential of the stem is as expected. We plot the water potential of the leaves along with it to see if they follow the same pattern.
plot(
plotgraph(sol, plantstructure, varname = :Ψ, structmod = [:Stem, :Leaf], title = "Low transpiration"),
plotgraph(sol2, plantstructure, varname = :Ψ, structmod = [:Stem, :Leaf], title = "High transpiration"),
ylims = (-0.35, -0.1), yaxis = "Ψ"
)
The plot shows us the expected pattern for both plant parts, verifying that the pre-built functionality can capture this basic hydraulic phenomenon. For more interesting applications and more advanced functionality, we refer to the subsequent tutorials.
The water potential behaving differently on the second day compared to the first is caused by changes in the metabolite concentration between those days, illustrated in the plots below.
plotgraph(sol, plantstructure, varname = :M, structmod = [:Leaf, :Stem])
plot(
plotgraph(sol, plantstructure, varname = :Π, structmod = :Leaf),
plotgraph(sol, plantstructure, varname = :P, structmod = :Leaf),
plotgraph(sol, plantstructure, varname = :Ψ, structmod = :Leaf,
xlabel = "Time (hr)"),
layout = (3, 1)
)