Tutorial 2: Package validation

In this tutorial, we cover some of the more advanced functionality from the package in order to recreate an hydraulic tree growth model from literature. The following topics are discussed:

Defining a new geometrical shape for the structural modules.
Defining a large plant structure with the use of graph rewriting (which can be considered a generalisation of L-systems).
Defining new functional modules.
- Using time-series data in a functional module.
- Defining a new, asymmetric connection module.
Using the advanced visualisation functionality to compare simulation outputs with data.

Setup

using Pkg; Pkg.activate("../..")

using PlantModules

using PlantGraphs, ModelingToolkit, OrdinaryDiffEq, Plots

using DataInterpolations, ForwardDiff

using PlutoUI; TableOfContents()

Context

In this tutorial, we recreate the dynamic sap flow model by Perämäki et al. (2001) in PlantModules.jl. The model simulates stem diameter variation in Scots pine trees caused by transpiration-driven changes in water tension and wood elasticity. The tree is modelled as a series of cylindrical stem and branch segments with a constant length and a radius that decreases as they move farther from the base of the tree. Water flow is modelled as a linear function of the water pressure difference between segments, including the effects of gravity, and the cross-sectional area connecting two segments. For boundary conditions, the outermost branch segments have a set transpiration rate per unit needle area based on measurements, and the base stem segment has a water inflow based on the water potential difference with the soil, of which the water potential is also based on measurements.

This model is an interesting case study for a few reasons:

It is also based on hydraulics-based plant growth, so we can use the hydraulics-based core functionality from the package.
The structure is composed of reasonably large number of segments from two different types, which lends itself well to a modular modeling approach.
There are plenty of interesting things happening both in the structural and functional definition that allow us to showcase the extensibility of the package.
The paper includes input- and output data, allowing for model validation.

Preparatory calculations

In order to replicate the model from the paper, we start by defining the structural and functional parameters used there. The paper was transparent enough that all functional parameters can simply be copied (with some unit conversions), though the structure of the tree was not provided exactly, so some estimations had to be made. The most influential of these on the model output is the fraction of the stem radius which consists of water-conducting sapwood (in contrast to the non-conducting heartwood), which was estimated visually.

begin
    # dimensions of compartments

    ## stem segments
    segment_length = 0.2 * 1e2 # cm; from m (input)
    tree_length = 13.2 * 1e2 # cm; from m (input)
    n_segments = tree_length ÷ segment_length

    ## stem radius at tree base
    stem_r0 = 0.2 * 1e2 / 2 # cm; from the diameter in m (input)
    radial_fraction_sapwood = 0.25 # (estimate)
    sapwood_r0 = stem_r0 * radial_fraction_sapwood

    ## crown
    crown_base_height = 6.0 * 1e2 # cm; from m (input)
    crown_length = tree_length - crown_base_height
    n_crown_segments = crown_length ÷ segment_length
    crown_start_age = (n_segments - n_crown_segments)

    ## branching
    branching_frequency = 3 # (estimate)
    n_branches = n_crown_segments ÷ branching_frequency

    Δr = sapwood_r0/(n_branches + 1) 
        # how much the radius of a stem segment decreases whenever the stem branches
        # we assume the tree decreases linearly in radius from the given starting radius at the bottom to a radius of near 0 at the top
        # (+1 to prevent stem radius from hitting 0)

    ## branch radius
    	# calculate using da vinci's rule: total cross area remains constant
        # before and after branching
    function branch_radius(r, Δr)
        area_0 = r^2*pi
        area_new = (r - Δr)^2*pi
        area_branch = area_0 - area_new
        return sqrt(area_branch / pi)
    end

    ## needle areas
    total_needle_area = 50.8 * 1e2^2 # cm^2; from m^2 (input)
    needle_areas = [i * total_needle_area / sum(1:n_branches)
                    for i in n_branches:-1:1]

    ## branchgrowth frequency
    branchgrowth_frequency = 3 
        # how many steps it takes branches to grow a new segment (estimate)

    ## root dimensions
    total_root_area = 96.1 * 1e2^2 # cm^2; from m^2 (input)

    # parameters

    ## hydraulic conductivities
    permeability = 2.5e-12 * 1e2^2 # cm^2; from m^2 (input)
    η = 1.0 * 1e-9 * (1/3600) # MPa hr; from mPa s (dynamic viscosity of water)
    l = segment_length # cm
    ρ_w = 1.0 # g / cm^3 (density of water)
    K_s_stem = ρ_w * permeability / η / l # g / hr / MPa / cm^2

    L_p = 3.2e-8 * 1e2 * 3600 # cm / hr / MPa; from m / s / MPa (input)
    K_roots = ρ_w * L_p * total_root_area # g / hr / MPa

    ## elastic modulus
    E_D_r_stem = 0.15 * 1e3 # MPa; from GPa (input)
    E_D_stem = [E_D_r_stem, 17.5 * E_D_r_stem] # Perämäki et al. use a constant ratio of 1/17.5 between the radial and longitudinal elastic moduli
end;

Structural definition

We introduce two new possibilities for the definition of plant structure.

Firstly, we need to define a new shape for our structural modules: the stem and branches consist of a hollow cylinder of growing, water-conducting sapwood with an "inert" core of heartwood in the middle. The base hydraulics module models a growing, water-conducting element, so we want to only model the sapwood.

Secondly, the model structure described in the paper consists of a reasonably large amount of segments (over a hundred) that adhere to some basic rules, so we define the structure by use of graph rewriting, using the functionality from PlantGraphs.jl.

Defining a new shape

In order to define a new shape, we define a new composite type as a subtype of the PlantModules.ModuleShape abstract type. Afterwards, we extend the following functions for the shape MyShape with dimensions the vector D, also listed in the PlantModules.ModuleShape docstring:

getdimensionality(m::MyShape): Return how many dimensions define the shape.
volume(m::MyShape, D): Calculate the volume of the shape.
cross_area(m::MyShape, D): Calculate the cross-sectional area of the shape.
surface_area(m::MyShape, D): Calculate the surface area of the shape.
- This one is used for the photosynthesis module, which we won't be using this tutorial, so we don't need to extend it.

"""
    HollowCylinder <: ModuleShape

A compartment shape representing a hollow cylinder. It is defined by two dimensions: the radius and the length, and the attribute `frac_sapwood`, denoting the fraction of the radius that is conducting sapwood.
"""
struct HollowCylinder <: PlantModules.ModuleShape
    frac_sapwood::Float64
end

PlantModules.getdimensionality(::HollowCylinder) = 2

PlantModules.cross_area(hc::HollowCylinder, D::AbstractArray) = D[1]^2 * π * (2 / hc.frac_sapwood - 1)

PlantModules.volume(hc::HollowCylinder, D::AbstractArray) = cross_area(hc, D) * D[2]

PlantModules.surface_area(hc::HollowCylinder, D::AbstractArray) = 2 * (D[1]/frac_sapwood) * π * D[2]

Plant structural module types

The plant structural modules are defined with an attribute for their dimensions and their height, as these are variables that need to differ between nodes. We define separate modules for the outermost segments of both the stem and the branches: this makes the graph rewriting easier, and the tips of the branches have a different functional behaviour from the rest of the branch, as all transpiration is assigned to them. The branch tips have a needle_area attribute which describes the total area of needles on that branch, which is used to calculate its transpiration.

struct Stem <: PlantGraphs.Node
    D::Vector
    h
end

struct Branch <: PlantGraphs.Node
    D::Vector
    h
end

struct StemTip <: PlantGraphs.Node
    D::Vector
    h
    branch_nr
    age
end

struct BranchTip <: PlantGraphs.Node
    D::Vector
    h
    age
    needle_area
end

Growing the plant

Graph rewriting is done by defining a set of rules that dictate how the graph needs to be rewritten every step and an initial graph. The rewriting rules are then applied a number of time to achieve the final graph.

PlantGraphs.jl implements their rules with three arguments:

The node type the rule applies to.
A function describing the condition for when to apply the rule.
A function describing how the node is transformed when the rule is applied.

This section is very context-specific and not functionality provided by this package, so we refer to the PlantGraphs.jl docs for specifics.

Rules

The outermost stem segment grows another thinner stem segment and a branch every branching_frequency steps, if the stem is old enough to form a crown:

branching_rule = Rule( 
    StemTip,
    lhs = st -> data(st).age >= crown_start_age && data(st).age % 
        branching_frequency == 0,
    rhs = st -> Stem(data(st).D, data(st).h) + (
        BranchTip([branch_radius(data(st).D[1], Δr), segment_length], 
                  data(st).h, 0, needle_areas[data(st).branch_nr]),
        StemTip(data(st).D - [Δr, 0.0], data(st).h + segment_length, 
                data(st).branch_nr + 1, data(st).age + 1)
    )
);

The outermost stem segment grows another stem segment if not branching:

vertical_growth_rule = Rule( 
    StemTip,
    lhs = st -> !(data(st).age >= crown_start_age && data(st).age % 
        branching_frequency == 0),
    rhs = st -> Stem(data(st).D, data(st).h) + 
        StemTip(data(st).D, data(st).h + segment_length, 
                data(st).branch_nr, data(st).age + 1)
);

The outermost branch segments grown another branch segment every branchgrowth_frequency steps:

branchgrowth_rule = Rule(
    BranchTip,
    lhs = br -> data(br).age % branchgrowth_frequency == 0,
    rhs = br -> Branch(data(br).D, data(br).h) +
        BranchTip(data(br).D, data(br).h, data(br).age + 1, 
                  data(br).needle_area)
);

The outermost branch segments also age if they are not growing:

branchage_rule = Rule(
    BranchTip,
    lhs = br -> !(data(br).age % branchgrowth_frequency == 0),
    rhs = br -> BranchTip(data(br).D, data(br).h, data(br).age + 1, 
                          data(br).needle_area)
);

Instantiation & rewriting

Finally, we define our initial graph and apply the rewriting rules. For the final step, we remove the StemTip node, as we don't require it.

axiom = StemTip([sapwood_r0, segment_length], segment_length/2, 1, 0);

begin
    plant = Graph(axiom = axiom, rules = 
        (vertical_growth_rule, branching_rule, branchage_rule, branchgrowth_rule)
    )
    for _ in 1:(n_segments)
        rewrite!(plant)
    end
    PlantGraphs.prune!(plant.graph, getid([node for node in getnodes(plant) if getstructmod(node) == :StemTip][1])) # remove StemTip
end

plotstructure(plant)

Connecting to the environment

We will again be using a single soil and air compartment, so there is little new in connecting the plant graph to the environment. However, for the visualisation of the resulting structure, removing the names of the structural modules with names = "" can give a clearer visualisation.

struct Soil <: PlantGraphs.Node end

struct Air <: PlantGraphs.Node end

graphs = [plant, Soil(), Air()];

intergraph_connections = [(1, 2) => (getnodes(plant)[1], :Soil), (1, 3) => (:BranchTip, :Air)];

plantstructure = PlantStructure(graphs, intergraph_connections);

plotstructure(plantstructure, names = "")

Function definition

New functional modules

As mentioned before, the model has a set transpiration from the branch tips to the air, defined per unit needle area, and a total needle area for every branch tip. To implement this, we need two new modules:

A module for the branch tips describing their total needle area.
A connection module between branch tips and the air that describes the transpiration in function of a given time-series for transpiration (per unit needle area) and the needle area of the branch tip segment.

Functionality is defined in terms of ModelingToolkit.jl Systems, so we recommend unfamiliar users to look at their docs, specifically the non-DSL way of defining systems.

We need our independent variable and derivative operator to be consistent with those used by the other functional modules, so we import then from the package first:

import PlantModules: t, d

Needle area module

The needle area of a branch is considered to be constant through time, so the functional module is pretty straightforward. The only tricky part is that we need to define it as a constant-valued variable rather than a parameter. This is because we need to reference this needle area in the transpiration connection, and we can only reference variables, not parameters.

function needle_area_module(; name, needle_area)
    @parameters begin
        needle_area_val = needle_area, 
            [description = "Value of total needle area on branch"]
    end 
    @variables begin
        needle_area(t), [description = "Total needle area on branch"]
    end
    eqs = [needle_area ~ needle_area_val]
    return System(eqs, t; name)
end;

Fixed transpiration

Next is the transpiration connection. We have time-series data for the transpiration per unit needle area, which we read in first.

readcsv(filepath) = readlines(filepath) .|>
    x -> split(x, ", ") .|> x -> parse(Float64, x);

transpiration_data = readcsv("./data/transpiration_data.csv");

begin
    transp_times = first.(transpiration_data)
    transp_values = last.(transpiration_data)
    scatter(transp_times, transp_values, xlabel = "Time (h)", ylabel = "Transpiration (mg / s / m^2)", legend = false)
end

To use the data in a functional module, we first need to correct the units to match the ones of our hydraulic module. Secondly, we need to turn it into a function that will return the transpiration at any timepoint, so some interpolation is required. We use a simple LinearInterpolation from the DataInterpolations.jl package for this.

correct_transpiration_units(x) = x * 1e-3 * 3600 * (1e-2)^2; # Unit conversion from (mg / s / m^2) to (g / hr / cm^2)

transpiration_rate = LinearInterpolation(
    correct_transpiration_units.(transp_values), 
    transp_times
);

The actual module definition of the fixed transpiration connection is rather complicated. We have two main complications: Firstly, this specific connection has no inherent directionality, so we need to define it manually. Secondly, we need to relate the variables used in the connection to the variables of the nodes it connects.

First, the specification of directionality. Let's compare with the usual hydraulic connection, for which water flow is defined in function of the water potential difference of the nodes it connects. Take the example of a connection between leaves and the air. If we consider the connection from the leaves' perspective or the air's perspective, we will both times have the same water potential difference but with a different sign. From both perspectives we will therefore also have the same water flow but with a different sign. This is perfect because our water flow from the air to the leaves needs to be the exact opposite of the flow from the leaves to the air. Because the flow direction is determined by the water potentials of the nodes it connects, we don't have to specify it ourselves.

Our fixed transpiration rate function, in contrast, does not give any indication what direction the water flow needs to go in, only the value. That is why we need to manually specify this when we define the connection. We can do this by using the special original_order keyword argument, which has the value true when the order between the structural modules of the node and the neighboring node is the same as originally specified by the user in the connecting_modules dictionary, and false otherwise. We'll here choose the order "branch first, air second". In this case, the water flow needs to be negative when the nodes are in the original order (water flows away from the tree during transpiration), and positive when not.

The second complication is inherent to all connection modules: we need to specify how the variables from our connection (i.e. the graph edge) relate to those of the structural modules it connects (i.e. the graph nodes). This is done by defining a function with the following inputs:

node_MTK: The functional module of the first node in the connection.
nb_node_MTK: The functional module of the second node in the connection.
connection_MTK: The functional module of the connection itself.
original_order (optional): Whether the first and second nodes are in the original order specified in the connecting_modules dictionary.

This function is then returned as the second output of our connection function. For our model, the variable that needs to be connected is the needle area. As it is only defined for the BranchTip nodes and not the Air node, we need to use the original_order argument to make sure the variable is connected to the system of the BranchTip node.

For completeness' sake: The water flux variable F also needs to be connected to the node systems. However, the way we usually want to connect water fluxes to the nodes is by defining the net water influx of a node as the sum of all connected water fluxes. Because the variable of a node is connected to the variables of multiple edge systems, this specification needs to be given separately. We don't need to change how water fluxes relate to the node systems for this model, but it is possible. See the PlantCoupling docstring for more details.

function fixed_transpiration_connection(; name, original_order)
    @variables begin
        F_s(t), [description = "Specific water flux from compartment 2 
            to compartment 1"]
        F(t), [description = "Water flux from compartment 2 to compartment 1"]
        needle_area(t), [description = "Total needle area on branch"]
    end

    polarity = original_order ? -1 : 1

    eqs = [
        F_s ~ transpiration_rate(t),
        F ~ polarity * F_s * needle_area,
    ]

    get_connection_eqset(node_MTK, nb_node_MTK, connection_MTK, original_order) = (
        original_order ? 
            [connection_MTK.needle_area ~ node_MTK.needle_area] : [connection_MTK.needle_area ~ nb_node_MTK.needle_area]
    )

    return System(eqs, t; name), get_connection_eqset
end;

Coupling

The coupling of function and structure remains, luckily, straightforward. We only have to pay special attention to the order of the BranchTip and Air in connecting_modules, as this will determine how the original_order variable works for our connection module.

module_coupling = Dict(
    :Stem => [hydraulic_module, constant_carbon_module, K_module],
    :Branch => [hydraulic_module, constant_carbon_module, K_module],
    :BranchTip => [hydraulic_module, needle_area_module,
                   constant_carbon_module, K_module],
    :Soil => [environmental_module, Ψ_soil_module, constant_K_module],
    :Air => [environmental_module],
);

connecting_modules = Dict(
    (:Soil, :Stem) => constant_hydraulic_connection,
    (:Stem, :Stem) => hydraulic_connection,
    (:Stem, :Branch) => hydraulic_connection,
    (:Branch, :Branch) => hydraulic_connection,
    (:Branch, :BranchTip) => hydraulic_connection,
    (:BranchTip, :Air) => fixed_transpiration_connection,
);

plantcoupling = PlantCoupling(; module_coupling, connecting_modules);

Parameters

The parameter specification mainly comes down to changing the default parameter values to those prescribed by the paper. Note that when we define new parameters or variables they need to be given a module-wide default value, even if those values are never used (here the case for needle_area). Some striking parameters are the zeros for the extensibility ϕ_D and metabolite concentration M: this was chosen because the model only considers elastic diameter changes and no irreversible growth, and does not consider carbon dynamics.

begin
    default_changes = Dict(
        :needle_area => 0.0,
        :E_D => E_D_stem, :K_s => K_s_stem,
        :ϕ_D => 0.0, :M => 0.0, 
        :Ψ => PlantModules.soilfunc(0.9), 
        :shape => HollowCylinder(radial_fraction_sapwood)
    )

    module_defaults = Dict(
        :Soil => Dict(:W_max => 1e6, :W_r => 0.9),
    )

    connection_values = Dict(
        (:Soil, :Stem) => Dict(:K => K_roots),
    )

    plantparams = PlantParameters(; default_changes, module_defaults,
                                  connection_values)
end;

Running the system

System generation and running also remains largely the same, with the exception that we specify to solve what solver we want. FBDF is a good solver for large, stiff systems of DAEs, as can be found on DifferentialEquations.jl's overview of ODE solvers.

system = generate_system(plantstructure, plantcoupling, plantparams);

tspan = (0.0, 24.0);

Warning

Problem construction of DAEs uses SciMLBase.SCCNonlinearProblem by default to solve the initialization problem. For the functional modules provided by PlantModules.jl, this currently takes much longer than the alternative option SciMLBase.NonlinearProblem for large systems. We can specify we want the alternative method for problem initialization using use_scc = false.

prob = ODEProblem(system, [], tspan, sparse = true, use_scc = false);

sol = solve(prob, FBDF());

Results

We now validate our model by comparing it to the simulation results and measured outputs found in the paper. Firstly, we verify that our transpiration was implemented correctly by plotting the total water influx of the air. After converting it to the units used in the paper, we can visually compare the results. Secondly, we compare our simulated water pressures for three different stem segments at different heights to those found in the paper. Finally, we compare our simulated diameter variation through time with the measured diameter variation from the paper.

Check transpiration

When we either want to plot very specific variables, or if we want to use ModelingToolkit.jl's plotting functionality to apply a transformation to variables before visualizing them, the get_subsystem_variables function can be used to extract the Symbolic representation of specific variables.

air_water_inflow = get_subsystem_variables(system, plantstructure, :ΣF, :Air)[1]

$$\begin{equation} \mathtt{Air158.{\Sigma}F}\left( t \right) \end{equation}$$

We can then define a function that changes the transpiration from the units used in our model (g / hr) to those in the paper (mg / m^2 / s).

transpiration_unit_conversion(var) = var * 1e3 / (total_needle_area * (1e-2)^2) / 3600;

plot(sol, idxs = [transpiration_unit_conversion(air_water_inflow)], label = false,
    xlabel = "Time of day (h)", ylabel = "Transpiration (mg / m^2 / s)",
    xticks = 0:3:24, ylims = (0, 10), yticks = 0:2:10, color = :black)

Water tension

We use the same approach to select the pressure variable from Stem nodes at specific heights of the tree.

begin
    stem_pressures = get_subsystem_variables(system, plantstructure, :P, :Stem)
    pressure_segment_heights = [10.0, crown_base_height, 1200.0]
    pressure_segment_nrs = ceil.(Int64, pressure_segment_heights / segment_length)
    
    plot(sol, idxs = stem_pressures[pressure_segment_nrs], ylims = (-0.4, 0.0),
         xticks = 0.0:3.0:24.0, xlabel = "Time of day (h)", color = :black,
         linewidth = [1 2 3], ylabel = "Water tension (MPa)",
         label = ["Base of stem" "Crown base" "Top of tree"])
end

Diameter variation

Finally, we can compare our simulation with actual measurements.

begin
    dimension_variables = get_subsystem_variables(system, plantstructure, :D, :Stem)

    diameter_change_mm(var, sol) = 2*10*(var - sol[var][1]) 
        # substract first value and change from radius in cm to diameter in mm
    diameter_segment_heights = [6.5 * 1e2, 2.5 * 1e2]
    diameter_segment_nrs = ceil.(Int64, diameter_segment_heights ./ segment_length)
    
    diameter_datas = [readcsv("./data/diameter_data_$(idx).csv") 
                      for idx in ("c", "d")]
    
    subplots = [plot(), plot()]
    for i in eachindex(subplots)
        plot!(subplots[i], sol, idxs = 
              [diameter_change_mm(
                  dimension_variables[diameter_segment_nrs[i]][1], sol
              )],
              label = "Simulated", color = :black, linewidth = 1)
        plot!(subplots[i], first.(diameter_datas[i]), last.(diameter_datas[i]), 
              label = "Data", color = :black, linewidth = 2)
        plot!(subplots[i], 
              title = "Stem at height $(diameter_segment_heights[i]) cm")
    end
    plot(subplots..., layout = (2, 1), ylabel = "Diameter change (mm)", 
         size = (800, 600), margins = 5*Plots.mm, ylims = (-0.07, 0.01), 
         yticks = -0.07:0.01:0.01, xticks = 0:3:24, xlabel = "Time of day (h)")
end

Uncertainty analysis

Perämäki et al. give a value range for the radial elasticity and the specific hydraulic conductivity of the tree segments. To inspect the effect of the uncertainty on these parameters on the results of our simulation, we conduct two types of uncertainty analysis here.

We will choose the diameter of the stem segment at 250 cm as our output variable of interest.

D1_at_250cm = dimension_variables[diameter_segment_nrs[2]][1];

Monte Carlo

First we perform a simple Monte Carlo experiment to visualize the effect of different values of the radial elasticity on the stem diameter variation.

We start by defining the radial elasticity's range of values and a function to sample from said range.

E_D_r_range = [0.03 * 1e3, 0.27 * 1e3] # MPa (from GPa)

2-element Vector{Float64}:
  30.0
 270.0

sample_range(a, b) = (b-a) * rand() + a

sample_range (generic function with 1 method)

Recreating the system and problem with generate_system and ODEProblem is a costly way to change the parameter values of the problem. We can circumvent this by using ModelingToolkit.jl's remake function to change the parameter values of the problem directly. PlantModules.jl's defines the helper function remake_graphsystem to make the remaking easier for systems with many subsystems.

begin
    p_montecarlo = plot()
                       
    for _ in 1:20
        E_D_r_sample = sample_range(E_D_r_range[1], E_D_r_range[2])
        E_D_sample = [E_D_r_sample, 17.5 * E_D_r_sample] 
            # Perämäki et al. use a constant ratio of 1/17.5 between the radial and longitudinal elastic moduli
    
        prob2 = remake_graphsystem(
            prob, system, plantstructure, :E_D, 
            [:Stem, :Branch, :BranchTip], E_D_sample
        )
        sol2 = solve(prob2, FBDF(), reltol = 1e-1) 
            # We use a higher relative tolerance for faster solving. Note that this can cause a SingularException error now and again but the faster solving time is generally worth it

        plot!(p_montecarlo, sol2, 
              idxs = [diameter_change_mm(D1_at_250cm, sol)], label = false, line_z = E_D_r_sample)
    end
    plot!(
        p_montecarlo, xlabel = "Time of day (h)", ylabel = "Diameter change (mm)",
        size = (800, 600), margins = 5*Plots.mm, ylims = (-0.12, 0.01), 
        yticks = -0.12:0.01:0.01, xticks = 0:3:24
    )
end

Local sensitivity

Using ForwardDiff.jl, it is possible to perform automatic differentiation to compute gradients through our simulation. We can apply this to investigate the effect of changes in the specific hydraulic conductivity on the stem diameter.

We begin by defining a function that returns our output variable of interest given the value of the specific hydraulic conductivity. remake_graphsystem only works when we change the value of a parameter to another one of the same type, which is not the case when we perform automatic differentiation. Instead we retrieve the specific hydraulic conductivities from all our stem and branch subsystems with the get_subsystem_variables function and then manually use pure remake to allow replacing parameters with values of different types, at the cost of performance.

K_s_vars = get_subsystem_variables(system, plantstructure, :K_s, [:Stem, :Branch]);

function get_diameter(new_K_s)
    newprob = remake(prob, p = Pair.(K_s_vars, [new_K_s]))
    newsol = solve(newprob, FBDF(), saveat = 0.1)
    return newsol[D1_at_250cm]
end

get_diameter (generic function with 1 method)

We can then call ForwardDiff.derivative on this function to get the derivative around a given value. We choose the value used in the simulation.

K_sensitivity = ForwardDiff.derivative(get_diameter, K_s_stem);

plot(
    tspan[1]:0.1:tspan[2], K_sensitivity, 
    xticks = 0:3:24, xlabel = "Time of day (h)", legend = false, 
    title = "Derivative of stem radius w.r.t. the hydraulic conductivity", 
    size = (800, 600), margins = 5*Plots.mm
)

The plot shows that an increase in hydraulic conductivity will have a positive impact on the diameter and that the magnitude of the effect follows the same pattern as the transpiration throughout the day, which intuitively makes sense!