Predator-prey palatability

This example compares two ways of configuring predator-prey palatability in the Agate.jl-NiPiZD model. By default, the model computes the palatability matrix from allometric traits based on the preferred predator:prey size ratio (optimum_predator_prey_ratio, abbreviated here as Vopt). We will change this default value to use a different optimum. Additionally, we will show how the allometric calculation can be bypassed by supplying an explicit custom palatability_matrix.

After inspecting the resulting interaction matrices, each configuration is run in a 0D box-model setup so their tracer trajectories can be compared.

Loading dependencies

The example uses Agate.jl, Oceananigans.jl, and OceanBioME.jl for the ecosystem simulation. CairoMakie.jl is used for plotting.

using Agate
using Agate.Introspection: interaction_matrix, tracer_names
using Agate.Library.Light
using OceanBioME
using OceanBioME: Biogeochemistry
using Oceananigans
using Oceananigans.Fields: FunctionField
using Oceananigans.Units
using CairoMakie

Construct the ecosystem configurations

default_bgc = Agate.Models.NiPiZD.construct()
default_pal = interaction_matrix(default_bgc, :palatability_matrix)

The default configuration uses the model's allometric palatability calculation. By default, predators prefer prey that are 10 times smaller than themselves (Vopt = 10).

Instead, we will change that optimum to Vopt = 5 for the two zooplankton consumers. Vector parameters can be overridden by plankton SizeClass name, so only the relevant entries need to be supplied. Omitted entries are filled from the model defaults.

vopt_bgc = Agate.Models.NiPiZD.construct(;
    parameters=(; optimum_predator_prey_ratio=(Z_1=5.0, Z_2=5.0))
)
vopt_pal = interaction_matrix(vopt_bgc, :palatability_matrix)

Another option is to provide the palatability matrix directly. This is useful for experiments where the intended interaction structure is known explicitly, rather than derived from size traits.

custom_bgc = Agate.Models.NiPiZD.construct(;
    palatability_matrix=[0.0 1.0; 1.0 0.0]
)
custom_pal = interaction_matrix(custom_bgc, :palatability_matrix)

Helper for labelled matrix plots

function plot_matrix!(fig, position, table; title)
    mat = Matrix(table.matrix)
    nrows, ncols = size(mat)

    ax = Axis(
        fig[position...];
        title,
        xlabel=string(table.column_axis),
        ylabel=string(table.row_axis),
        xticks=(1:ncols, string.(table.columns)),
        yticks=(1:nrows, string.(table.rows)),
        yreversed=true,
        aspect=DataAspect(),
    )

    hm = heatmap!(ax, 1:ncols, 1:nrows, mat'; colorrange=(0, 1))

    for row in 1:nrows, col in 1:ncols
        text!(
            ax,
            col,
            row;
            text=string(round(mat[row, col]; digits=2)),
            align=(:center, :center),
            fontsize=12,
        )
    end

    return hm
end

Compare the palatability matrices

matrix_fig = Figure(; size=(1100, 380), fontsize=16)

hm_pal_default = plot_matrix!(matrix_fig, (1, 1), default_pal; title="Vopt = 10 (default)")

plot_matrix!(matrix_fig, (1, 2), vopt_pal; title="Vopt = 5")

plot_matrix!(matrix_fig, (1, 3), custom_pal; title="Custom")

Colorbar(matrix_fig[1, 4], hm_pal_default; label="palatability")

Label(matrix_fig[0, 1:4], "Predator-prey palatability"; fontsize=22)

save("predator_prey_palatability_matrices.png", matrix_fig; px_per_unit=1)

matrix_fig

Run each configuration in a box model

function run_box_model(bgc, filename)
    grid = BoxModelGrid()
    clock = Clock(; time=zero(grid))
    PAR = FunctionField{Center,Center,Center}(CyclicalPAR(-10), grid; clock)
    light_attenuation = PrescribedPhotosyntheticallyActiveRadiation(PAR)
    bgc_model = Biogeochemistry(bgc; light_attenuation)
    full_model = BoxModel(; grid, clock, biogeochemistry=bgc_model)

    set!(full_model; N=7.0, P_1=0.01, Z_1=0.01, P_2=0.1, Z_2=0.01, D=0.01)

    simulation = Simulation(full_model; Δt=240minutes, stop_time=1095days)

    simulation.output_writers[:fields] = JLD2Writer(
        full_model,
        full_model.fields;
        filename,
        schedule=TimeInterval(1day),
        overwrite_existing=true,
    )

    run!(simulation)

    return filename
end

default_filename = run_box_model(default_bgc, "predator_prey_palatability_default.jld2")
vopt_filename = run_box_model(vopt_bgc, "predator_prey_palatability_vopt5.jld2")
custom_filename = run_box_model(custom_bgc, "predator_prey_palatability_custom.jld2")

[ Info: Initializing simulation...
┌ Warning: error ArgumentError("a group or dataset named grid is already present within this group") thrown when trying to serialize the grid at serialized/grid
└ @ Oceananigans.OutputWriters ~/.julia/packages/Oceananigans/9oDU5/src/OutputWriters/jld2_writer.jl:251
[ Info:     ... simulation initialization complete (1.453 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (555.880 ms).
[ Info: Simulation is stopping after running for 5.640 seconds.
[ Info: Simulation time 1095 days equals or exceeds stop time 1095 days.
[ Info: Initializing simulation...
┌ Warning: error ArgumentError("a group or dataset named grid is already present within this group") thrown when trying to serialize the grid at serialized/grid
└ @ Oceananigans.OutputWriters ~/.julia/packages/Oceananigans/9oDU5/src/OutputWriters/jld2_writer.jl:251
[ Info:     ... simulation initialization complete (10.958 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (102.282 μs).
[ Info: Simulation is stopping after running for 3.605 seconds.
[ Info: Simulation time 1095 days equals or exceeds stop time 1095 days.
[ Info: Initializing simulation...
┌ Warning: error ArgumentError("a group or dataset named grid is already present within this group") thrown when trying to serialize the grid at serialized/grid
└ @ Oceananigans.OutputWriters ~/.julia/packages/Oceananigans/9oDU5/src/OutputWriters/jld2_writer.jl:251
[ Info:     ... simulation initialization complete (11.263 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (108.012 μs).
[ Info: Simulation is stopping after running for 3.665 seconds.
[ Info: Simulation time 1095 days equals or exceeds stop time 1095 days.

Load the output time series

tracer_syms = tracer_names(default_bgc)

default_timeseries = (;
    (s => FieldTimeSeries(default_filename, string(s))[1, 1, 1, :] for s in tracer_syms)...
)

vopt_timeseries = (;
    (s => FieldTimeSeries(vopt_filename, string(s))[1, 1, 1, :] for s in tracer_syms)...
)

custom_timeseries = (;
    (s => FieldTimeSeries(custom_filename, string(s))[1, 1, 1, :] for s in tracer_syms)...
)

times = FieldTimeSeries(default_filename, string(first(tracer_syms))).times

Compare the box-model output

simulation_fig = Figure(; size=(900, 900), fontsize=16)

Label(simulation_fig[0, 1:2], "Box-model response"; fontsize=22)

for (idx, sym) in enumerate(tracer_syms)
    row = 1 + div(idx - 1, 2)
    col = 1 + mod(idx - 1, 2)

    ax = Axis(
        simulation_fig[row, col];
        ylabel=string(sym),
        xlabel="Days",
        title="$(sym) concentration (mmol N / m³)",
    )

    lines!(
        ax,
        times / day,
        getproperty(default_timeseries, sym);
        label="Vopt = 10",
        linewidth=2,
    )
    lines!(
        ax,
        times / day,
        getproperty(vopt_timeseries, sym);
        label="Vopt = 5",
        linewidth=2,
    )
    lines!(
        ax,
        times / day,
        getproperty(custom_timeseries, sym);
        label="custom",
        linewidth=2,
    )
    axislegend(ax; position=:rt)
end

save("predator_prey_palatability_simulation.png", simulation_fig; px_per_unit=1)

simulation_fig

This page was generated using Literate.jl.