Mutualism simulation

How to run

There are three things you can do with this code:

• Run a Single simulation
• Run a competition parameter sweep
• Generate plots

Before the code will run, a few file directory names must be updated. In tauSweep.m, change outputDir to the directory where MATLAB should save all matfiles. The script will create this directory if it does not exist.

In plotjob.sh, the first argument given to the generatePlots function should be the directory where the sweep’s matfiles were saved. The second argument should be the directory where MATLAB should save the figures generated by the function.

Running single simulations

To run a single simulation with specific parameter values, simply cd into the project directory and call the function simMutualism(). Add any custom parameter values as arguments to the function, with the variable name in quotes, followed by the desired value, e.g.:

simMutualism('tau12', 0.2, 'tau21', 0.3, 'sigma_sq', 0.1)

The specified parameters are used in the filename for the resulting mat file that the function generates.

ODE parameters

The possible variables to modify for the ODE parameters are rP, rF1, rF2, alphaPF1, alphaPF2, alphaF1P, alphaF2P, qP, qF1, qF2, betaP, betaF1, betaF2, cP, cF1, cF2, dP, dF1, dF2, hPF1, hPF2, hF1P, hF2P, eP, eF1, eF2, deltaP, deltaF1, deltaF2, tau12, tau21.

Running competition parameter sweeps

To run a competition parameter sweep, use the function tauSweep().

Generating plots

The generatePlots() function will create a set of plots from a competition parameter sweep. To generate individual plots for a single simulation, use plotFinalPopSpace(), plotPopSpaceTime(), plotSpeedTime(), and plotRangeTime().

Using git to keep track of parameter sweep settings

One issue with running several sweeps with slightly different parameter values is that it’s easy to lose track of the parameter values you used for a given sweep. This can cause problems if you need to reproduce a certain set of results, or re-run a sweep with slightly different parameter values. It’s possible to use git branches to keep track of all the different sweeps we run.

We have our main git branch, which can be considered the default settings and parameter values. If we need to change some parameter values for a specific sweep, we likely don’t want to save these changes to our main branch since it would no longer represent the default settings. When we want to run a new sweep, we can use git to create and checkout a new branch:

git checkout -b SWEEP-BRANCH-NAME

The command git checkout allows you to switch between your branches, and in this case with the -b option, create a new one and switch to it.

The idea here is that we name our branch something meaningful that tells us what parameter values we used for our sweep. For example, if I’m running a sweep with the useDeltaDispKernels parameter set to true, I could name the branch delta-disp-kernels.

Then, I’d go into the tauSweep.m file and change the outputDir variable to ~/delta-disp-kernels so my results are saved into a new directory. I’d also add the necessary arguments to the simMutualism function call (e.g. setting useDeltaDispKernels to true). I would commit these changes, which only saves them to the new delta-disp-kernels branch.

After running the simulation, I could switch back to the main branch with the command

git checkout main

The delta-disp-kernels branch is saved if I ever need to check what parameter values I used for the sweep or if I need to re-run it. But since I haven’t modified the main branch at all, the default settings are preserved. I can now create another branch if I need to run a different sweep with different parameter values.

How the code works

Simulation

Here we simulate the outcome of many generations of growth-dispersal cycles for three mutualist species.

function simMutualism(varargin)

% we save the given parameters to a variable so we can use them to
% name files and label plots
parameters = varargin;

Parameters

Simulation parameters

The number of growth-dispersal cycles (i.e. iterations) will change based on whether a steady state is reached for all three species. A species reaches a steady state when the variance between the last 10 spread speed values is at or below some threshold.

So the number of iterations defined here is just the minimum number of iterations that a simulation will run. We also set the maximum number of iterations, so even if a steady state isn’t reached, the simulation will end after this many cycles. This value is also used to allocate space in a number of arrays used in the simulation. It’s faster (although more costly in memory) to preallocate the maximum number of rows we may need to use, rather than to try and increase the size of several large arrays.

We also define the number of additional iterations (iterationStep) to add to the simulation before checking for steady states again.

%% for simulation

p = inputParser;
p.KeepUnmatched = true;

% minimum number of cycles of growth and dispersal
addParameter(p, 'iterations', 50, @isnumeric);
addParameter(p, 'maxIterations', 480, @isnumeric);
addParameter(p, 'iterationStep', 50, @isnumeric);
addParameter(p, 'outputDir', './', @isfolder);
addParameter(p, 'steadyStateThreshold', 1e-03, @isnumeric);
addParameter(p, 'diameter', 1200, @isnumeric);
addParameter(p, 'sigma_sq', 0.25, @isnumeric); % Dispersal variance
addParameter(p, 'deltaP', 0.0, @isnumeric);
addParameter(p, 'deltaF1', 0.9, @isnumeric);
addParameter(p, 'deltaF2', 0.1, @isnumeric);
addParameter(p, 'useDeltaDispKernels', false, @islogical);


parse(p, varargin{:});

% I wish I knew a better way to get rid of all the p.Results that get attached
% inputParser parameters
iterations = p.Results.iterations;
maxIterations = p.Results.maxIterations;
sigma_sq = p.Results.sigma_sq;
deltaP = p.Results.deltaP;
deltaF1 = p.Results.deltaF1;
deltaF2 = p.Results.deltaF2;

if iterations > maxIterations
    disp("Warning: the value of iterations is greater than or ");
    disp("equal to maxIterations, so maxIterations has been increased.");
    maxIterations = iterations;
end

iterationStep = p.Results.iterationStep;
outputDir = p.Results.outputDir;
steadyStateThreshold = p.Results.steadyStateThreshold;

%total size of landscape along positive x-axis (so half the total landscape)
diameter = p.Results.diameter;

Space parameters

Here we create the one-dimensional landscape in which the species will disperse.

linspace(x1, x2, n) creates a vector of n points between points x1 and x2. Spacing between points is (x2-x1)/(n-1). See linspace documentation.

%% Initialize space parameters
lowval = 1e-9;
nodes = (2^16) + 1; %total points in space -- 65537
radius = diameter / 2;
x = linspace(-radius, radius, nodes);
x2 = linspace(-diameter, diameter, 2 * nodes - 1);
dx = diameter / (nodes - 1);

Initialization

Here we initialize the arrays we’ll use throughout the simulation. The speed arrays save the instantaneous or average spread speed of a species for each iteration. The range edge arrays keep track of the furthest spatial location on one side of a species’ range. The n population arrays keep track of each species’ population density across the entire spatial range.

% preallocate arrays for max possible iterations + 1
[instantSpeedP, avgSpeedP, instantSpeedF1, avgSpeedF1, instantSpeedF2, avgSpeedF2] = deal(zeros(1, maxIterations + 1));

[rangeEdgeP,rangeEdgeF1, rangeEdgeF2] = deal(zeros(1, maxIterations + 1));

[nP, nF1, nF2] = deal(zeros(maxIterations + 1, length(x)));

Dispersal kernels

We use a Gaussian dispersal kernel for each species. At some point we’ll rewrite this to allow for other dispersal kernel functions.

If we want to have the dependence parameter affect the dispersal kernel, we can use the useDeltaDispKernels with the value true when calling simMutualism(). Otherwise, the default is to have dispersal unaffected by dependence.

if p.Results.useDeltaDispKernels
    % gaussian dispersal kernels
    kP = exp(-(x2 .^ 2) / (2 * sigma_sq)) ./ sqrt(2 * pi * sigma_sq);
    kF1 = exp(-(x2 .^ 2) / (2 * sigma_sq * deltaF1)) ./ sqrt(2 * pi * sigma_sq * deltaF1);
    kF2 = exp(-(x2 .^ 2) / (2 * sigma_sq * deltaF2)) ./ sqrt(2 * pi * sigma_sq * deltaF2);
else
    kP = exp(-(x2 .^ 2) / (2 * sigma_sq)) ./ sqrt(2 * pi * sigma_sq);
    kF1 = exp(-(x2 .^ 2) / (2 * sigma_sq)) ./ sqrt(2 * pi * sigma_sq);
    kF2 = exp(-(x2 .^ 2) / (2 * sigma_sq)) ./ sqrt(2 * pi * sigma_sq);
end

Initial population densities

We set the initial population densities across the spatial range.

% SET THE INITIAL CONDITIONS
irad = 2; % Initial condition range
initDensities = [0.1,0.1,0.1];
nThreshold = 0.05; % critical threshold for edge of wave
temp_P = find(abs(x) <= irad); %locate all values in the array x that lie b/w +irad and -irad units of space
temp_F1 = find(abs(x) <= irad);
temp_F2 = find(abs(x) <= irad);

nP(1,temp_P) = initDensities(1) * normpdf(x(temp_P),0,1); %Computes pdf values evaluated at the values in x i.e. all x(temp) values for the normal distribution with mean 0 and standard deviation 1.
nF1(1,temp_F1) = initDensities(2) * normpdf(x(temp_F1),0,1);
nF2(1,temp_F2) = initDensities(3) * normpdf(x(temp_F2),0,1);

Initial front location

% FIND THE INITIAL FRONT LOCATION
% find the farthest distance travelled by the population above a certain threshold density and assign it to front
frontP = find(nP(1,:) >= nThreshold,1,'last');
frontF1 = find(nF1(1,:) >= nThreshold,1,'last');
frontF2 = find(nF2(1,:) >= nThreshold,1,'last');

% the initial front is obtained from initialization which will be in the first
% row of 'n'
if frontP
  rangeEdgeP(1) = interp1(nP(1,frontP:frontP+1),x(frontP:frontP+1),nThreshold);
end
if frontF1
  rangeEdgeF1(1) = interp1(nF1(1,frontF1:frontF1+1),x(frontF1:frontF1+1),nThreshold);
end

if frontF2
  rangeEdgeF2(1) = interp1(nF2(1,frontF2:frontF2+1),x(frontF2:frontF2+1),nThreshold);
end

Simulating growth and dispersal over many generations

generation = 1;
%% Looping for growth and dispersal
while generation <= iterations

Growth phase

% for ode45
tspan = [0, 10];

%Growth
y0 = [nP(generation,:);nF1(generation,:);nF2(generation,:)];

% reshape happens such that 3 consecutive rows for nP, nF1, and nF2 values
% are stacked
y0 = reshape(y0, 3*length(y0), 1);

%remember to alter where the dep_p and dep_f are being called from
[t,y] = ode45(@(t,y) growthODEs(t,y, varargin{:}), tspan, y0);


% We just want the results of the growth phase (end)
fP = y(end,(1:3:end)); % final row; element 1, +3, elem. 4, etc. until end
fF1 = y(end,(2:3:end));
fF2 = y(end,(3:3:end));

Dispersal phase

%   DISPERSAL
    n1P = fft_conv(kP,fP);
    n1F1 = fft_conv(kF1,fF1);
    n1F2 = fft_conv(kF2,fF2);

    nP(generation + 1,:) = dx*n1P(nodes:length(x2));
    nF1(generation + 1,:) = dx*n1F1(nodes:length(x2));
    nF2(generation + 1,:) = dx*n1F2(nodes:length(x2));

    nP(generation + 1,1) = nP(generation + 1,1)/2;
    nP(generation + 1,nodes) = nP(generation + 1,nodes)/2;

    nF1(generation + 1,1) = nF1(generation + 1,1)/2;
    nF1(generation + 1,nodes) = nF1(generation + 1,nodes)/2;

    nF2(generation + 1,1) = nF2(generation + 1,1)/2;
    nF2(generation + 1,nodes) = nF2(generation + 1,nodes)/2;

    % gives location of random places where numbers are above zero due to some
    % numerical errors
    temp_P = find(nP(generation + 1,:) < lowval);
    temp_F1 = find(nF1(generation + 1,:) < lowval);
    temp_F2 = find(nF2(generation + 1,:) < lowval);

    % set the places with those numerical errors to zero
    nP(generation + 1,temp_P) = zeros(size(nP(generation + 1,temp_P)));
    nF1(generation + 1,temp_F1) = zeros(size(nF1(generation + 1,temp_F1)));
    nF2(generation + 1,temp_F2) = zeros(size(nF2(generation + 1,temp_F2)));

    frontP = find(nP(generation + 1,:) >= nThreshold,1,'last');
    frontF1 = find(nF1(generation + 1,:) >= nThreshold,1,'last');
    frontF2 = find(nF2(generation + 1,:) >= nThreshold,1,'last');

    % if any of the species' range edge is equal to the edge of the entire
    % spatial range, stop the growth-dispersal loop. We set total iterations to
    % the last iteration + 1 so the data is still usable.
    if (frontP == nodes) | (frontF1 == nodes) | (frontF2 == nodes)
        error("Warning: the simulation has stopped because the edge of the landscape was reached.");
    end

    if frontP
         rangeEdgeP(generation + 1) = interp1(nP(generation + 1,frontP:frontP + 1),x(frontP:frontP + 1), nThreshold);
    end

    if frontF1
         rangeEdgeF1(generation + 1) = interp1(nF1(generation + 1, frontF1:frontF1 + 1), x(frontF1:frontF1 + 1), nThreshold);
    end

    if frontF2
         rangeEdgeF2(generation + 1) = interp1(nF2(generation + 1,frontF2:frontF2 + 1), x(frontF2:frontF2 + 1), nThreshold);
    end

    %latest position of wave edge - initial position of wave edge divided by time
    avgSpeedP(generation) = (rangeEdgeP(generation + 1) - rangeEdgeP(1)) / generation;

    instantSpeedP(generation) = rangeEdgeP(generation + 1) - rangeEdgeP(generation);

    instantSpeedF1(generation) = rangeEdgeF1(generation + 1) - rangeEdgeF1(generation);

    %latest position of wave edge - initial position of wave edge divided by time
    avgSpeedF1(generation) = (rangeEdgeF1(generation + 1) - rangeEdgeF1(1)) / generation;

    %latest position of wave edge - initial position of wave edge divided by time
    instantSpeedF2(generation) = rangeEdgeF2(generation + 1) - rangeEdgeF2(generation);
    avgSpeedF2(generation) = (rangeEdgeF2(generation + 1) - rangeEdgeF2(1)) / generation;

Determine whether to continue running the simulation for more iterations

    % check for steady state, and determine whether to run for more generations
    if (generation == iterations)

        % if not all species at steady state
        if ~(isSpeciesSteadyState(instantSpeedP, steadyStateThreshold, generation) && isSpeciesSteadyState(instantSpeedF1, steadyStateThreshold, generation) && isSpeciesSteadyState(instantSpeedF2, steadyStateThreshold, generation))

            % end the simulation if you've hit maxIterations
            if generation == maxIterations
                error("Warning: The simulation for tau12 = %s and tau21 = %s has reached the maxIterations value of %s.", p.Results.tau12, p.Results.tau21, maxIterations)
            end

            % iterations close to the max
            if iterations >= (maxIterations - iterationStep)
                iterations = maxIterations;
            else
                iterations = iterations + iterationStep;
            end
        end
    end

    generation = generation + 1;

% while loop end
end

Checking if a species is at a steady state

This function takes the spread speed values for a given species and checks to see if the variance in the last 10 values is at or below a threshold to determine whether a steady state has been reached.

function isSteadyState = isSpeciesSteadyState(speed, tolerance, generation)
% takes a matrix of speed values and checks whether the variance in the last 10 values is at or below a threshold

    variance = sqrt(var(speed((generation - 9):generation)));

    if variance <= tolerance
        isSteadyState = true;
    else
        isSteadyState = false;
    end
end

Which simulations never reached a steady state?

Sometimes a simulation ends because the maxIterations number is reached, rather than actually reaching a steady state. Here we write a function to tell us if a simulation ended because it reached the maxIterations. It will iterate through all the files in the given directory and print a list of the parameter regimes for which a steady state wasn’t reached.

function getNoSteadyState(sweepDir)

    files = dir(fullfile(sweepDir, '*.mat'));

    for file = 1:length(files)
        curFile = matfile(fullfile(sweepDir, files(file).name));

        parameters = curFile.parameters;

        % get the values of tau12 and tau21
        tau12 = parameters{find(strcmp('tau12', parameters)) + 1};
        tau21 = parameters{find(strcmp('tau21', parameters)) + 1};

        if curFile.iterations == curFile.maxIterations
            disp(strcat("The simulation of tau12 = ", num2str(tau12, "%.2f"), " and tau21 = ", num2str(tau21, "%.2f"), " reached the maxIterations value of ", num2str(curFile.maxIterations)));
        else
            disp(strcat("The simulation of tau12 = ", num2str(tau12, "%.2f"), " and tau21 = ", num2str(tau21, "%.2f"), " ran for ", num2str(curFile.iterations), " iterations"))
        end

        clear curFile;
    end
end

Generate and save a mat file for the simulation

Using maxIterations to create the initial arrays means that these arrays may be storing many more rows than is actually necessary. Since we’re saving these to mat files, we can reduce the size before saving by resizing the arrays. By getting rid of extra rows, we can also use the end index to get the population densities of the final iteration.

Then we can save our results to a mat file, which can then be used to generate figures, identify outcomes, etc. The filename string can be reused for saving figures as well. It takes any explicitly defined parameters from the call to simMutualism() and appends the names and values to filename.

%% Save a mat file with the current parameter values

nP = nP(1:(iterations + 1), :);
nF1 = nF1(1:(iterations + 1), :);
nF2 = nF2(1:(iterations + 1), :);

instantSpeedP = instantSpeedP(1, 1:(iterations + 1));
instantSpeedF1 = instantSpeedF1(1, 1:(iterations + 1));
instantSpeedF2 = instantSpeedF2(1, 1:(iterations + 1));

% classify outcome here so we don't have to do it later
outcome = classifyOutcome(nF1(end,:), nF2(end,:), nThreshold);

%% Save a mat file with the current parameter values

filename = 'results';
formatSpec = '%.2f';

if ~(isempty(parameters))
    for i = 1:length(parameters)
        param = parameters{i};

        if isnumeric(param)
            param = num2str(param, formatSpec);
        elseif strcmp(param, 'outputDir') || islogical(param) || isfolder(param)
            continue
        else
            param = string(param);
        end

        filename = strcat(filename, '_', param);
    end
end

filename = strcat(filename, '.mat');

save(strcat(outputDir, filename), 'nP', 'nF1', 'nF2', 'iterations', 'nThreshold', 'instantSpeedP', 'instantSpeedF1', 'instantSpeedF2', 'filename', 'parameters', 'x', 'maxIterations', 'diameter', 'outcome');

% end of simMutualism function
end

Growth equations function

Here we define the growth of each species using a system of ODEs.

System of Equations (growthODEs.m)

Function definition

With varargin, we can optionally use parameter values other than the defaults, e.g. growthODEs(t, y, 'rP', 0.4). We need to use an inputParser to manage the function’s parameters.

function dydt = growthODEs(t, y, varargin)

Default parameter values

We set our default parameter values here. If the parameter is not explicitly defined in the function call, then these default values are used.

%% Default ODE parameter values

default_nodes = (2^16) + 1;

% intrinsic growth
default_rP = 0.3;
default_rF1 = 0.3;
default_rF2 = 0.3;

% mutualism benefits
default_alphaPF1 = 0.01;
default_alphaPF2 = 0.01;
default_alphaF1P = 0.5;
default_alphaF2P = 0.5;

default_qP = 1.0;
default_qF1 = 1.0;
default_qF2 = 1.0;

% mutualism costs
default_betaP = 0.0;
default_betaF1 = 0.0;
default_betaF2 = 0.0;

default_cP = 1.0;
default_cF1 = 1.0;
default_cF2 = 1.0;

% death rate
default_dP = 0.1;
default_dF1 = 0.1;
default_dF2 = 0.1;

% saturation
default_hPF1 = 0.3;
default_hPF2 = 0.3;
default_hF1P = 0.3;
default_hF2P = 0.3;

default_eP = 0.3;
default_eF1 = 0.3;
default_eF2 = 0.3;

% = 0.0;
default_deltaP = 0.1;
default_deltaF1 = 0.9;
default_deltaF2 = 0.1;

% competition: tau12 is the effect F2 has on F1; tau21 is effect of F1 on F2
default_tau12 = 0.0;
default_tau21 = 0.0;

Adding parameters with inputParser

See inputParser and addParameter documentation. By setting p.KeepUnmatched = true, we can pass along all the parameters given in the simMutualism function call and just ignore the ones that are not relevant to the ODE parameters.

p = inputParser;
p.KeepUnmatched = true;

addRequired(p, 't');
addRequired(p, 'y');

%% Optional ODE parameters

addParameter(p, 'nodes', default_nodes);

% intrinsic growth rates
addParameter(p, 'rP', default_rP);
addParameter(p, 'rF1', default_rF1);
addParameter(p, 'rF2', default_rF2);

% mutualism benefits
addParameter(p, 'alphaPF1', default_alphaPF1);
addParameter(p, 'alphaPF2', default_alphaPF2);
addParameter(p, 'alphaF1P', default_alphaF1P);
addParameter(p, 'alphaF2P', default_alphaF2P);

addParameter(p, 'qP', default_qP );
addParameter(p, 'qF1', default_qF1);
addParameter(p, 'qF2', default_qF2);

% mutualism costs
addParameter(p, 'betaP', default_betaP);
addParameter(p, 'betaF1', default_betaF1);
addParameter(p, 'betaF2', default_betaF2);

addParameter(p, 'cP', default_cP);
addParameter(p, 'cF1', default_cF1);
addParameter(p, 'cF2', default_cF2);

% death rate
addParameter(p, 'dP', default_dP);
addParameter(p, 'dF1', default_dF1);
addParameter(p, 'dF2', default_dF2);

% saturation
addParameter(p, 'hPF1', default_hPF1);
addParameter(p, 'hPF2', default_hPF2);
addParameter(p, 'hF1P', default_hF1P);
addParameter(p, 'hF2P', default_hF2P);

addParameter(p, 'eP', default_eP);
addParameter(p, 'eF1', default_eF1);
addParameter(p, 'eF2', default_eF2);

% mutualism dependence
addParameter(p, 'deltaP', default_deltaP);
addParameter(p, 'deltaF1', default_deltaF1);
addParameter(p, 'deltaF2', default_deltaF2);

% competition
addParameter(p, 'tau12', default_tau12);
addParameter(p, 'tau21', default_tau21);

parse(p, t, y, varargin{:});

% relabel variables so they're easier to read in the equation

t = p.Results.t;
y = p.Results.y;
nodes = p.Results.nodes;

% intrinsic growth
rP = p.Results.rP;
rF1 = p.Results.rF1;
rF2 = p.Results.rF2;

% mutualism benefits
alphaPF1 = p.Results.alphaPF1;
alphaPF2 = p.Results.alphaPF2;
alphaF1P = p.Results.alphaF1P;
alphaF2P = p.Results.alphaF2P;

cP = p.Results.cP;
cF1 = p.Results.cF1;
cF2 = p.Results.cF2;

% death rate
dP = p.Results.dP;
dF1 = p.Results.dF1;
dF2 = p.Results.dF2;

% saturation
hPF1 = p.Results.hPF1;
hPF2 = p.Results.hPF2;
hF1P = p.Results.hF1P;
hF2P = p.Results.hF2P;

% mutualism dependence
deltaF1 = p.Results.deltaF1;
deltaF2 = p.Results.deltaF2;

% competition: tau12 is the effect F2 has on F1; tau21 is effect of F1 on F2
tau12 = p.Results.tau12;
tau21 = p.Results.tau21;

y = reshape(y,3,nodes);
dydt  = zeros(size(y));

Species P

$$\begin{align*} \frac{dP}{dt} = P [ r_P + \left( c_1 \left[\frac{\alpha_{PF_1} F_1}{h_{PF_1} + F_1} + \frac{\alpha_{PF_2} F_2}{h_{PF_2} + F_2} \right] \right) - d_P P ] \end{align*}$$

% rename variables so equations are easier to read
P = y(1,:);
F1 = y(2,:);
F2 = y(3,:);

dydt(1,:) = P .* (rP + (cP * (alphaPF1 .* F1 ./ (hPF1 + F1) + alphaPF2 .* F2 ./ (hPF2 + F2))) - dP .* P);

Species F₁

$$\begin{align*} \frac{dF_1}{dt} = F_1[(1 - \delta_{F_1})r_{F_1} + \delta_{F_1} \left( c_2 \left[\frac{\alpha_{F_1} P}{h_{F_1} + P} \right] \right) - \tau_{12} F_2 - d_{F_1} F_1] \end{align*}$$

dydt(2,:) = F1 .* ((1 - deltaF1) * rF1 + deltaF1 * (cF1 * (alphaF1P .* P) ./ (hF1P + P)) - (tau12 .* F2) - dF1 .* F1);

Species F₂

$$\begin{align*} \frac{dF_2}{dt} = F_2[(1 - \delta_{F_2}) r_{F_2} + \delta_{F_2} \left(c_2 \left[\frac{\alpha_{F_2} P}{h_{F_2} + P} \right] \right) - \tau_{21} F_1 - d_{F_2} F_2] \end{align*}$$

dydt(3,:) = F2 .* ((1 - deltaF2) * rF2 + deltaF2 * (cF2 * (alphaF2P .* P) ./ (hF2P + P)) - (tau21 .* F1) - dF2 .* F2);

Reshape

    dydt = reshape(dydt,3*nodes,1);
end

Parameter sweep

Sweep script

This is the main file to be edited when running parameter sweeps. The simMutualism() function requires an output directory as an argument. It can take any ODE parameter as an optional argument. to override a default value, use the parameter variable name then a value, i.e. simMutualism(outputDir, 'tau12', 0.3, 'tau21', 0.14). These variables will get added to the filename of the exported mat file at the end of the simulation.

We can adjust values of maxIterations in the for loop to allow for longer simulations of $τ$ values that we know will take longer to reach a steady state. By allow for more iterations for only these values, we limit the number of very large mat files.

% use integers for the number of iterations to run (rather than the actual
% values of tau12 and tau21) because it seems parfor requires it

rangeStep = 0.01;

outputDir = '/home/shawa/lutzx119/deltaDispSweep/';

mkdir(outputDir)

% instead of using a for loop for the tau12 values, we can use Slurm to set up
% jobs for each tau12 value. To change the range of tau12 values, modify the
% "SBATCH --array=" line in the Slurm job script.
tau12 = rangeStep * str2num(getenv("SLURM_ARRAY_TASK_ID"));

parfor j = 0:40

    tau21 = j * rangeStep;
    simMutualism('outputDir', outputDir, 'tau12', tau12, 'tau21', tau21, 'useDeltaDispKernels', true);
end

Slurm job script

The SBATCH lines must be at the top of the script. Anything before that will cause an error with Slurm.

#!/bin/bash -l
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=4
#SBATCH --mem-per-cpu=4G
#SBATCH --time=4:00:00
#SBATCH --array=0-40
#SBATCH --mail-type=ALL
#SBATCH [email protected]
#SBATCH --output=/home/shawa/lutzx119/reports/tausweep-%j.out

cd /home/shawa/lutzx119/mutualism || return
module purge

module load matlab
matlab -nodisplay <tauSweep.m

Using the Slurm --array command

Note the #SBATCH --array=0-40 command. For a parameter sweep, we might normally use two for-loops to iterate through a range of values for two parameters (in this case $τ_ { 12 }$ and $τ_ { 21 }$). Instead, we can replace the outer for loop with Slurm’s --array command. It takes a range of numbers (or a comma separated list in brackets, like [5, 10, 15, 25]), and creates a separate, parallelized task for each one. We access the task ID with the environment variable SLURM_ARRAY_TASK_ID.

In our case, we want to run simulations for all values of $τ_ { 12 }$ and $τ_ { 21 }$ in the range 0 - 0.4 with a step size of 0.01. Since the Slurm array command only recognizes integers, we use --array=0-40. Then in our Matlab sweep script, we remove our outer for-loop and replace wherever we were using the for-loop index variable with str2num(getenv("SLURM_ARRAY_TASK_ID")) * 0.01. Now Slurm will set up one job with 41 sub-tasks that run in parallel, one for each value of $τ_ { 12 }$.

Making the right job request

This introductory guide from Princeton Research Computing was very helpful. They also have specific instructions for Matlab.

The key takeaway is that requesting more resources as a way to speed up a job is usually a bad idea. Often it’ll get stuck in the queue and any performance gains are offset by this wait. Even worse, Matlab doesn’t typically benefit from multiple nodes/ntasks. It’s best to simply use --nodes=1 and --ntasks=1. Using the --array command as mentioned above with automatically spread the parameter sweep simulations across many CPUs/cores, so there’s no need to do anything else except request sufficient memory, either with --mem= or --mem-per-cpu=.

It is possible to parallelize your code as well, but it’s not totally clear whether this is always beneficial. You can replace a for-loop with parfor, and then add --cpus-per-task to your Slurm script. This is another case where more is not always better, --cpus-per-task=4 has given me the best results in my very informal testing. More CPUs and your job ends up in the queue for much longer. Now for each array sub-task, 4 CPUs will divide up the inner for-loop in the sweep script. The problem is that these extra CPUs can mean that your other array sub-tasks get stuck in the queue since you’re using more resources per task.

Function to classify outcome (classifyOutcome.m)

This function takes the final population densities of species $F_1$ and $F_2$ and classifies the outcome of the simulation. The possible outcomes are:

• F1 dominance (outcome = 1)
• F2 dominance (outcome = 2)
• local coexistence (outcome = 3)
• local coexistence with F1 dominance (outcome = 4)
• local coexistence with F2 dominance (outcome = 5)
• regional coexistence (outcome = 6)
• unknown (outcome = 7).

First we find the values above nThreshold across the landscape — this gives us each species final range. We use the max function to determine whether $F_1$ or $F_2$ had the bigger range.

In order to make classification easier we create a variable that tells us whether or not $F_1$ had the larger range than $F_2$, based on the result of the max function.

%% Function to classify outcome of a given simulation
function outcome = classifyOutcome(finalNF1, finalNF2, nThreshold)

    % get the ranges where F1 and F2 populations are above the threshold
    rangeF1 = find(finalNF1 >= nThreshold);
    rangeF2 = find(finalNF2 >= nThreshold);

    lenMaxRange = max(length(rangeF1), length(rangeF2));

    % range where one species exists but not the other
    exclusiveRange = setxor(rangeF1, rangeF2);

In the simplest cases, there were no population values above nThreshold for either $F_1$ or $F_2$; this means the other species competitively excluded it and we can classify the outcome as $F_1$ or $F_2$ dominance.

% if F2 is below the threshold across the total range, then classify as
% F1 dominance
if isempty(rangeF2)
    outcome = 1; % F1 dominance

% if F1 is below the threshold across the total range, then classify as
% F2 dominance
elseif isempty(rangeF1)
    outcome = 2; % F2 dominance

Next, we look to see if the $F$ species with the larger range was dominant for less than 0.05 of its total range. The setxor function gives us the areas of space where the species with the greater range competitively excluded the other. We determine the total length of these areas and then divide by maxRange to get the proportion of the total range where this species was dominant. If this proportion is less than the (arbitrary) threshold of 0.05, we classify this as local coexistence.

% find the range of values in rangeF1 or rangeF2 but not both
% if the proportion of this range over the total range is less than
% the arbitrary value 0.05, we call it local coexistence
elseif length(exclusiveRange)/lenMaxRange < 0.05
    outcome = 3; % Local coexistence

It’s possible that the proportion of space where the dominant species competitively excluded the other is greater than 0.05. In this case, we first determine if $F_1$ was the dominant species (i.e. it had the larger range). Since we’ve already found outcomes where the lengths of the ranges of $F_1$ and $F_2$ differ by less than 5%, we know that any outcomes found here will have at least some local dominance.

Since setxor(rangeF1, rangeF2) gives us any area of the landscape where one species competitively excluded the other, we use intersect to see if any of those areas fall within rangeF2. In other words, if $F_2$ competitively excluded $F_1$ for any proportion of the landscape. If so, we classify this as regional coexistence.

If not, this means that $F_1$ has regions of its total range where it has competitively excluded $F_2$ (the proportion of which must be greater than or equal to 0.05). We know from the comparisons above, however, that $F_1$ still occupies some proportion of the landscape, so we classify this as local coexistence with $F_1$ dominance.

We then make the same comparisons when $F_2$ has the larger range. Finally, we classify any outcome that does not fall into these categories as “unknown”, which most likely indicates some sort of error.

    elseif length(rangeF1) > length(rangeF2)

        % no F2 dominance
        if isempty(intersect(rangeF2, exclusiveRange))
            outcome = 4; % Local coexistence + F1 dominance
        % we find at least some F2 dominance
        else
            outcome = 6; % regional coexistence
        end

    elseif length(rangeF2) > length(rangeF1)

        % no F1 dominance
        if isempty(intersect(rangeF1, exclusiveRange))
            outcome = 5; % Local coexistence + F2 dominance
        else
            % we find at least some F1 dominance
            outcome = 6; % regional coexistence
        end
    else
        outcome = 7; % unknown
    end
end

Figures

Generate figures from paper

The parameter space plot always includes all the parameter values in the sweep. For the other plots, the variables tau12Range and tau21Range define for what range of parameter values the plots are generated.

function generatePlots(sweepDir, figDir, varargin)

    defaultTau12Range = 0.13:0.01:0.31;
    defaultTau21Range = 0.0:0.01:0.4;

    p = inputParser;
    addRequired(p, 'sweepDir', @isfolder);
    addRequired(p, 'figDir');
    addParameter(p, 'plotOutcomes', false, @islogical);
    addParameter(p, 'plotPopSpaceTime', false, @islogical);
    addParameter(p, 'plotFinalPopSpace', false, @islogical);
    addParameter(p, 'plotSpeedTime', false, @islogical);
    addParameter(p, 'tau12Range', defaultTau12Range, @isvector);
    addParameter(p, 'tau21Range', defaultTau21Range, @isvector);
    addParameter(p, 'taus', [], @ismatrix);

    parse(p, sweepDir, figDir, varargin{:});

    mkdir(figDir)

    if p.Results.plotOutcomes
        % get the heatmap of all the outcomes
        disp('Generating outcomes plot...')
        if isfolder(figDir)
            plotOutcomes(sweepDir, 'figDir', figDir);
        else
            error("figDir is not a folder")
        end
    end

    if p.Results.plotPopSpaceTime || p.Results.plotFinalPopSpace || p.Results.plotSpeedTime

        tau12Range = p.Results.tau12Range;
        tau21Range = p.Results.tau21Range;
        taus = p.Results.taus;

        % check to make sure generatePlots is given either tau ranges or pairs but not both
        if ~(isequal(tau12Range, defaultTau12Range) && isequal(tau21Range, defaultTau21Range)) && ~isempty(taus)

            error("Specify values for tau ranges or a vector of tau pair values, but not both")
        end

        if isempty(taus)
            for tau12 = tau12Range

                taus = [taus; ones(numel(tau21Range), 1) * tau12, tau21Range(:)];

            end
        end

        for i = 1:length(taus)

            formatSpec = '%.2f';

            % probably a better way to do this with regexp
            targetFile = dir(fullfile(sweepDir, strcat("*tau12_", num2str(taus(i, 1), formatSpec), "*tau21_", num2str(taus(i, 2), formatSpec), "*.mat")));

            filename = fullfile(sweepDir, targetFile.name);

            curFile = load(filename, 'iterations', 'filename', 'nP', 'nF1', 'nF2', 'nThreshold', 'x', 'instantSpeedP', 'instantSpeedF1', 'instantSpeedF2');

            if p.Results.plotPopSpaceTime
                plotPopSpaceTime(curFile, 'figDir', figDir);
            end

            if p.Results.plotFinalPopSpace
                plotFinalPopSpace(curFile, 'figDir', figDir);
            end

            if p.Results.plotSpeedTime
                plotSpeedTime(curFile, 'figDir', figDir);
            end

            clear curFile;
        end
    end
end

3D population density vs. space vs. time plot

These plots are helpful to see how the population densities change over time, but the 2D final spatial outcome plots are a little easier to read if all we care about is what happens at the steady state.

We generate a plot for each species, and they’re superimposed in a single figure.

function plotPopSpaceTime(simMatFile, varargin)

    p = inputParser;
    addRequired(p, 'simMatFile');
    addOptional(p,'createFile', true, @islogical);
    addOptional(p, 'figDir', './', @isfolder);

    parse(p, simMatFile, varargin{:});

    filename = simMatFile.filename;
    iterations = simMatFile.iterations;
    nP = simMatFile.nP;
    nF1 = simMatFile.nF1;
    nF2 = simMatFile.nF2;
    diameter = simMatFile.diameter;
    nThreshold = simMatFile.nThreshold;
    x = simMatFile.x;

    timeStep = round(iterations / 15);

    %% Figure for species P

    % if you're creating a file, don't display the figure in a window
    if p.Results.createFile
        f = figure('visible', 'off');
    else
        figure(1);
    end

    [xx,tt] = meshgrid(x,0:iterations);
    nlow = nP;
    nlow(nP >= nThreshold) = NaN;
    nP(nP < nThreshold) = NaN;

    rangeP = x(find(nP(end,:) >= nThreshold));

    rangeMin = min(rangeP);
    rangeMax = max(rangeP);

    hold on
    for i = 1:timeStep:iterations
        lineP = plot3(xx(i,:),tt(i,:),nP(i,:),'b', 'LineWidth', 3.0);
        plot3(xx(i,:),tt(i,:),nlow(i,:),'Color',0.8*[1 1 1]);
        grid on
    end
    % plot3(rangeEdgeP(1:11),0:10,nThreshold*ones(1,11),'k');
    axis([(rangeMin - 5) (rangeMax + 5) 0 iterations 0 6.25]);
    xticks([rangeMin 0 rangeMax]);
    xticklabels({num2str(-diameter/2), '0', num2str(diameter/2)})
    xlabel('Spatial range');
    ylabel('Generations');
    zlabel('Population density');
    % title('Species P');
    view(30,30);

    %% Figure for species F1
    [xx,tt] = meshgrid(x,0:iterations);
    nlow = nF1;
    nlow(nF1 >= nThreshold) = NaN;
    nF1(nF1 < nThreshold) = NaN;
    hold on
    for i = 2:timeStep:iterations
        lineF1 = plot3(xx(i,:),tt(i,:),nF1(i,:),'r','LineWidth', 3.0);
        plot3(xx(i,:),tt(i,:),nlow(i,:),'Color',0.8*[1 1 1]);
        grid on
    end

    %% Figure for species F2
    [xx,tt] = meshgrid(x,0:iterations);
    nlow = nF2;
    nlow(nF2 >= nThreshold) = NaN;
    nF2(nF2 < nThreshold) = NaN;
    hold on
    for i = 3:timeStep:iterations
        lineF2 = plot3(xx(i,:),tt(i,:),nF2(i,:),'g', 'LineWidth', 3.0);
        plot3(xx(i,:),tt(i,:),nlow(i,:),'Color',0.8*[1 1 1]);
        grid on
    end
    hold off

    legend([lineP lineF1 lineF2], {'P', 'F_1', 'F_2'});

    if p.Results.createFile
        [~, filename, ~] = fileparts(filename);
        filename = strcat('pop_space_time_', filename);
        savefig(strcat(p.Results.figDir, filename, '.fig'));
        saveas(strcat(p.Results.figDir, filename, '.png'));
        clf;
    end
end

Speed vs. time plot

function plotSpeedTime(simMatFile, varargin)

    p = inputParser;
    addRequired(p, 'simMatFile');
    addOptional(p,'createFile', true, @islogical);
    addOptional(p, 'figDir', './', @isfolder);

    parse(p, simMatFile, varargin{:});

    filename = simMatFile.filename;
    iterations = simMatFile.iterations;
    instantSpeedP = simMatFile.instantSpeedP;
    instantSpeedF1 = simMatFile.instantSpeedF1;
    instantSpeedF2 = simMatFile.instantSpeedF2;

    if p.Results.createFile
        f = figure('visible', 'off');
    else
        figure(1);
    end

    plot(1:iterations, instantSpeedP(1:iterations), 1:iterations, instantSpeedF1(1:iterations), 1:iterations, instantSpeedF2(1:iterations));
    legend('P', 'F1', 'F2');
    title(strcat(['Spread speed vs. time']));
    xlabel('iterations');
    ylabel('speed');

    if p.Results.createFile
        [~, filename, ~] = fileparts(filename);
        filename = fullfile(p.Results.figDir, strcat('speed_time_', filename));
        saveas(f, strcat(filename, '.fig'));
        saveas(f, strcat(filename, '.png'));
    end
end

Final population densities across space plot

%
% PLOTFINALPOPSPACE Plot the species' final spatial ranges from a single simulation.
%   PLOTFINALPOPSPACE(simMatFile) takes data from a matfile loaded with the matfile() function and creates a fig file and a png file in the current directory.
%
%   PLOTFINALPOPSPACE(simMatFile, 'figDir', './someDirectory/') sets the directory where the fig and png files are saved.
%
%   PLOTFINALPOPSPACE(simMatFile, 'createFile', false) does not save any files but instead displays the plot in a new window.
%
%   See also PLOTPOPSPACETIME, PLOTSPEEDTIME, PLOTRANGETIME.
function plotFinalPopSpace(simMatFile, varargin)

    p = inputParser;
    addRequired(p, 'simMatFile');
    addOptional(p,'createFile', true, @islogical);
    addOptional(p, 'figDir', './', @isfolder);

    parse(p, simMatFile, varargin{:});

    diameter = simMatFile.diameter;
    nP = simMatFile.nP;
    nF1 = simMatFile.nF1;
    nF2 = simMatFile.nF2;
    nThreshold = simMatFile.nThreshold;
    x = simMatFile.x;

    iterations = simMatFile.iterations;
    filename = simMatFile.filename;

    if p.Results.createFile
        f = figure('visible', 'off');
    else
        figure(1);
    end

    rangeP = find(nP(iterations,:) >= nThreshold);

    rangeMin = min(rangeP);
    rangeMax = max(rangeP);

    f.Position = [1 1 996 996];
    axis square;

    hold on
    plot(nP(iterations,:), LineWidth=1.5);
    plot(nF1(iterations,:), LineWidth=1.5);
    plot(nF2(iterations,:), LineWidth=1.5);
    xlim([(rangeMin - 1000) (rangeMax + 1000)]);
    % xticks([(rangeMin - 1000) ((rangeMin - 1000) * 2) (rangeMax + 1000)]);
    xlabel('Spatial range');
    ylabel('Population density');
    xticks([(rangeMin - 1000) (width(nP)/2) (rangeMax + 1000)]);
    xticklabels({num2str(int32(diameter*2/width(nP)*(rangeMin - 1000) - diameter)), '0', num2str(int32(diameter*2/width(nP)*(rangeMax + 1000) - diameter))});

    legend('P', 'F1', 'F2');
    hold off

    if p.Results.createFile
        [~, filename, ~] = fileparts(filename);
        filename = fullfile(p.Results.figDir, strcat('final_pop_space_', filename));
        saveas(f, strcat(filename, '.fig'));
        saveas(f, strcat(filename, '.png'));
        clf;
    end
end

Sweep outcomes plot

This function generates a heatmap of the outcomes of a $τ$ parameter sweep (it might be possible to make this more generic for other types of parameter sweeps in the future). It requires a directory where it can find mat files (the results of each simulation).

It can optionally take arguments to specify the range of values used in the parameter sweep (by default it assumes that we used the range 0:0.01:0.40 for both $τ_ { 12 }$ and $τ_ { 21 }$).

function plotOutcomes(sweepDir, varargin)

    p = inputParser;

    addRequired(p, 'sweepDir', @isfolder);
    addParameter(p, 'tau12Range', 0:0.01:0.40);
    addParameter(p, 'tau21Range', 0:0.01:0.40);
    addParameter(p, 'figDir', './', @isfolder);
    parse(p, sweepDir, varargin{:});

    tau12Range = p.Results.tau12Range;
    tau21Range = p.Results.tau21Range;
    figDir = p.Results.figDir;

    outcomes = zeros(length(tau12Range), length(tau21Range));

    files = dir(fullfile(sweepDir, '*.mat'));

    for file = 1:length(files)

        curFile = matfile(fullfile(sweepDir, files(file).name));

        parameters = curFile.parameters;
        % get the values of tau12 and tau21
        tau12 = parameters{find(strcmp('tau12', parameters)) + 1};
        tau21 = parameters{find(strcmp('tau21', parameters)) + 1};

        disp(strcat("The outcome of tau12 = ", num2str(tau12, "%.2f"), " and tau21 = ", num2str(tau21, "%.2f"), " is ", num2str(curFile.outcome)));

        % You can't use == for comparison of floating point numbers, you have to
        % use this ismembertol function The default tolerance is fine for this
        % purpose.
        outcomes(ismembertol(tau12Range, tau12), ismembertol(tau21Range, tau21)) = curFile.outcome;

        clear curFile;

    end

    f = figure('visible', 'off');
    heatmap(tau12Range, fliplr(tau21Range), rot90(outcomes));
    xlabel('tau_{12}');
    ylabel('tau_{21}');

    filename = fullfile(figDir, 'tauSweepOutcomesPlot');
    disp("Saving outcomes plot to %s...", filename)
    saveas(f, strcat(filename, '.fig'));
    saveas(f, strcat(filename, '.png'));

end