```
using VPL
using Base.Threads: @threads
using Plots
import Random
using FastGaussQuadrature
using Distributions
Random.seed!(123456789)
```

`Random.TaskLocalRNG()`

In this example we extend the forest example to have more complex, time- depedent development and growth based on carbon allocation. For simplicity, the model assumes a constant relative growth rate at the plant level to compute the biomass increment. In the next example this assumption is relaxed by a model of radiation use efficiency. When modelling growth from carbon allocation, the biomass of each organ is then translated in to an area or volume and the dimensions of the organs are updated accordingly (assuming a particular shape).

The following packages are needed:

```
using VPL
using Base.Threads: @threads
using Plots
import Random
using FastGaussQuadrature
using Distributions
Random.seed!(123456789)
```

`Random.TaskLocalRNG()`

The data types needed to simulate the trees are given in the following module. The differences with respect to the previous example are:

```
- Meristems do not produce phytomers every day
- A relative sink strength approach is used to allocate biomass to organs
- The geometry of the organs is updated based on the new biomass
- Bud break probability is a function of distance to apical meristem
```

```
# Data types
module TreeTypes
using VPL
using Distributions
# Meristem
Base.@kwdef mutable struct Meristem <: VPL.Node
::Int64 = 0 # Age of the meristem
ageend
# Bud
struct Bud <: VPL.Node end
# Node
struct Node <: VPL.Node end
# BudNode
struct BudNode <: VPL.Node end
# Internode (needs to be mutable to allow for changes over time)
Base.@kwdef mutable struct Internode <: VPL.Node
::Int64 = 0 # Age of the internode
age::Float64 = 0.0 # Initial biomass
biomass::Float64 = 0.0 # Internodes
length::Float64 = 0.0 # Internodes
width::Exponential{Float64} = Exponential(5)
sinkend
# Leaf
Base.@kwdef mutable struct Leaf <: VPL.Node
::Int64 = 0 # Age of the leaf
age::Float64 = 0.0 # Initial biomass
biomass::Float64 = 0.0 # Leaves
length::Float64 = 0.0 # Leaves
width::Beta{Float64} = Beta(2,5)
sinkend
# Graph-level variables -> mutable because we need to modify them during growth
Base.@kwdef mutable struct treeparams
# Variables
::Float64 = 2e-3 # Current total biomass (g)
biomass# Parameters
::Float64 = 1.0 # Relative growth rate (1/d)
RGR::Float64 = 1e-3 # Initial biomass of an internode (g)
IB0::Float64 = 0.1e6 # Specific internode weight (g/m3)
SIW::Float64 = 15.0 # Internode shape parameter (length/width)
IS::Float64 = 1e-3 # Initial biomass of a leaf
LB0::Float64 = 100.0 # Specific leaf weight (g/m2)
SLW::Float64 = 3.0 # Leaf shape parameter (length/width)
LS::Float64 = 1/0.5 # Bud break probability coefficient (in 1/m)
budbreak::Int64 = 5 # Number of days between phytomer production
plastochron::Float64 = 15.0 # Number of days that a leaf expands
leaf_expansion::Float64 = 140.0
phyllotaxis::Float64 = 30.0
leaf_angle::Float64 = 45.0
branch_angleend
end
import .TreeTypes
```

The methods for creating the geometry and color of the tree are the same as in the previous example.

```
# Create geometry + color for the internodes
function VPL.feed!(turtle::Turtle, i::TreeTypes.Internode, vars)
# Rotate turtle around the head to implement elliptical phyllotaxis
rh!(turtle, vars.phyllotaxis)
HollowCylinder!(turtle, length = i.length, height = i.width, width = i.width,
= true, color = RGB(0.5,0.4,0.0))
move return nothing
end
# Create geometry + color for the leaves
function VPL.feed!(turtle::Turtle, l::TreeTypes.Leaf, vars)
# Rotate turtle around the arm for insertion angle
ra!(turtle, -vars.leaf_angle)
# Generate the leaf
Ellipse!(turtle, length = l.length, width = l.width, move = false,
= RGB(0.2,0.6,0.2))
color # Rotate turtle back to original direction
ra!(turtle, vars.leaf_angle)
return nothing
end
# Insertion angle for the bud nodes
function VPL.feed!(turtle::Turtle, b::TreeTypes.BudNode, vars)
# Rotate turtle around the arm for insertion angle
ra!(turtle, -vars.branch_angle)
end
```

The meristem rule is now parameterized by the initial states of the leaves and internodes and will only be triggered every X days where X is the plastochron.

```
# Create right side of the growth rule (parameterized by the initial states
# of the leaves and internodes)
function create_meristem_rule(vleaf, vint)
= Rule(TreeTypes.Meristem,
meristem_rule = mer -> mod(data(mer).age, vars(mer).plastochron) == 0,
lhs = mer -> TreeTypes.Node() +
rhs Bud(),
(TreeTypes.Leaf(biomass = vleaf.biomass,
TreeTypes.= vleaf.length,
length = vleaf.width)) +
width Internode(biomass = vint.biomass,
TreeTypes.= vint.length,
length = vint.width) +
width Meristem())
TreeTypes.end
```

`create_meristem_rule (generic function with 1 method)`

The bud break probability is now a function of distance to the apical meristem rather than the number of internodes. An adhoc traversal is used to compute this length of the main branch a bud belongs to (ignoring the lateral branches).

```
# Compute the probability that a bud breaks as function of distance to the meristem
function prob_break(bud)
# We move to parent node in the branch where the bud was created
= parent(bud)
node # Extract the first internode
= filter(x -> data(x) isa TreeTypes.Internode, children(node))[1]
child = data(child)
data_child # We measure the length of the branch until we find the meristem
= 0.0
distance while !isa(data_child, TreeTypes.Meristem)
# If we encounter an internode, store the length and move to the next node
if data_child isa TreeTypes.Internode
+= data_child.length
distance = children(child)[1]
child = data(child)
data_child # If we encounter a node, extract the next internode
elseif data_child isa TreeTypes.Node
= filter(x -> data(x) isa TreeTypes.Internode, children(child))[1]
child = data(child)
data_child else
error("Should be Internode, Node or Meristem")
end
end
# Compute the probability of bud break as function of distance and
# make stochastic decision
= min(1.0, distance*vars(bud).budbreak)
prob return rand() < prob
end
# Branch rule parameterized by initial states of internodes
function create_branch_rule(vint)
= Rule(TreeTypes.Bud,
branch_rule = prob_break,
lhs = bud -> TreeTypes.BudNode() +
rhs Internode(biomass = vint.biomass,
TreeTypes.= vint.length,
length = vint.width) +
width Meristem())
TreeTypes.end
```

`create_branch_rule (generic function with 1 method)`

We need some functions to compute the length and width of a leaf or internode from its biomass

```
function leaf_dims(biomass, vars)
= biomass
leaf_biomass = biomass/vars.SLW
leaf_area = sqrt(leaf_area*4*vars.LS/pi)
leaf_length = leaf_length/vars.LS
leaf_width return leaf_length, leaf_width
end
function int_dims(biomass, vars)
= biomass
int_biomass = biomass/vars.SIW
int_volume = cbrt(int_volume*4*vars.IS^2/pi)
int_length = int_length/vars.IS
int_width return int_length, int_width
end
```

`int_dims (generic function with 1 method)`

Each day, the total biomass of the tree is updated using a simple RGR formula and the increment of biomass is distributed across the organs proportionally to their relative sink strength (of leaves or internodes).

The sink strength of leaves is modelled with a beta distribution scaled to the `leaf_expansion`

argument that determines the duration of leaf growth, whereas for the internodes it follows a negative exponential distribution. The `pdf`

function computes the probability density of each distribution which is taken as proportional to the sink strength (the model is actually source-limited since we imposed a particular growth rate).

```
sink_strength(leaf, vars) = leaf.age > vars.leaf_expansion ? 0.0 :
pdf(leaf.sink, leaf.age/vars.leaf_expansion)/100.0
plot(0:1:50, x -> sink_strength(TreeTypes.Leaf(age = x), TreeTypes.treeparams()),
= "Age", ylabel = "Sink strength", label = "Leaf") xlabel
```

```
sink_strength(int) = pdf(int.sink, int.age)
plot!(0:1:50, x -> sink_strength(TreeTypes.Internode(age = x)), label = "Internode")
```

Now we need a function that updates the biomass of the tree, allocates it to the different organs and updates the dimensions of said organs. For simplicity, we create the functions `leaves()`

and `internodes()`

that will apply the queries to the tree required to extract said nodes:

```
get_leaves(tree) = apply(tree, Query(TreeTypes.Leaf))
get_internodes(tree) = apply(tree, Query(TreeTypes.Internode))
```

`get_internodes (generic function with 1 method)`

The age of the different organs is updated every time step:

```
function age!(all_leaves, all_internodes, all_meristems)
for leaf in all_leaves
+= 1
leaf.age end
for int in all_internodes
+= 1
int.age end
for mer in all_meristems
+= 1
mer.age end
return nothing
end
```

`age! (generic function with 1 method)`

The daily growth is allocated to different organs proportional to their sink strength.

```
function grow!(tree, all_leaves, all_internodes)
# Compute total biomass increment
= vars(tree)
tvars = tvars.RGR*tvars.biomass
ΔB += ΔB
tvars.biomass # Total sink strength
= 0.0
total_sink for leaf in all_leaves
+= sink_strength(leaf, tvars)
total_sink end
for int in all_internodes
+= sink_strength(int)
total_sink end
# Allocate biomass to leaves and internodes
for leaf in all_leaves
+= ΔB*sink_strength(leaf, tvars)/total_sink
leaf.biomass end
for int in all_internodes
+= ΔB*sink_strength(int)/total_sink
int.biomass end
return nothing
end
```

`grow! (generic function with 1 method)`

Finally, we need to update the dimensions of the organs. The leaf dimensions are

```
function size_leaves!(all_leaves, tvars)
for leaf in all_leaves
= leaf_dims(leaf.biomass, tvars)
leaf.length, leaf.width end
return nothing
end
function size_internodes!(all_internodes, tvars)
for int in all_internodes
= int_dims(int.biomass, tvars)
int.length, int.width end
return nothing
end
```

`size_internodes! (generic function with 1 method)`

All the growth and developmental functions are combined together into a daily step function that updates the forest by iterating over the different trees in parallel.

```
get_meristems(tree) = apply(tree, Query(TreeTypes.Meristem))
function daily_step!(forest)
@threads for tree in forest
# Retrieve all the relevant organs
= get_leaves(tree)
all_leaves = get_internodes(tree)
all_internodes = get_meristems(tree)
all_meristems # Update the age of the organs
age!(all_leaves, all_internodes, all_meristems)
# Grow the tree
grow!(tree, all_leaves, all_internodes)
= vars(tree)
tvars size_leaves!(all_leaves, tvars)
size_internodes!(all_internodes, tvars)
# Developmental rules
rewrite!(tree)
end
end
```

`daily_step! (generic function with 1 method)`

The trees are initialized in a regular grid with random values for the initial orientation and RGR:

```
= rand(Normal(0.3,0.01), 10, 10)
RGRs histogram(vec(RGRs))
```

```
= [rand()*360.0 for i = 1:2.0:20.0, j = 1:2.0:20.0]
orientations histogram(vec(orientations))
```

`= [Vec(i,j,0) for i = 1:2.0:20.0, j = 1:2.0:20.0]; origins `

The following initalizes a tree based on the origin, orientation and RGR:

```
function create_tree(origin, orientation, RGR)
# Initial state and parameters of the tree
= TreeTypes.treeparams(RGR = RGR)
vars # Initial states of the leaves
= leaf_dims(vars.LB0, vars)
leaf_length, leaf_width = (biomass = vars.LB0, length = leaf_length, width = leaf_width)
vleaf # Initial states of the internodes
= int_dims(vars.LB0, vars)
int_length, int_width = (biomass = vars.IB0, length = int_length, width = int_width)
vint # Growth rules
= create_meristem_rule(vleaf, vint)
meristem_rule = create_branch_rule(vint)
branch_rule = T(origin) + RH(orientation) +
axiom Internode(biomass = vint.biomass,
TreeTypes.= vint.length,
length = vint.width) +
width Meristem()
TreeTypes.= Graph(axiom = axiom, rules = (meristem_rule, branch_rule),
tree = vars)
vars return tree
end
```

`create_tree (generic function with 1 method)`

As in the previous example, it makes sense to visualize the forest with a soil tile beneath it. Unlike in the previous example, we will construct the soil tile using a dedicated graph and generate a `Scene`

object which can later be merged with the rest of scene generated in daily step:

```
Base.@kwdef struct Soil <: VPL.Node
::Float64
length::Float64
widthend
function VPL.feed!(turtle::Turtle, s::Soil, vars)
Rectangle!(turtle, length = s.length, width = s.width, color = RGB(255/255, 236/255, 179/255))
end
= RA(-90.0) + T(Vec(0.0, 10.0, 0.0)) + # Moves into position
soil_graph Soil(length = 20.0, width = 20.0) # Draws the soil tile
= Scene(Graph(axiom = soil_graph));
soil render(soil, axes = false)
```

And the following function renders the entire scene (notice that we need to use `display()`

to force the rendering of the scene when called within a loop or a function):

```
function render_forest(forest, soil)
= Scene(vec(forest)) # create scene from forest
scene = Scene([scene, soil]) # merges the two scenes
scene render(scene)
end
```

`render_forest (generic function with 1 method)`

We may want to extract some information at the canopy level such as LAI. This is best achieved with a query:

```
function get_LAI(forest)
= 0.0
LAI @threads for tree in forest
for leaf in get_leaves(tree)
+= leaf.length*leaf.width*pi/4.0
LAI end
end
return LAI/400.0
end
```

`get_LAI (generic function with 1 method)`

We can now create a forest of trees on a regular grid:

```
= create_tree.(origins, orientations, RGRs);
forest render_forest(forest, soil)
for i in 1:50
daily_step!(forest)
end
render_forest(forest, soil)
```

And compute the leaf area index:

`get_LAI(forest)`

`0.24217455687553385`