% Supplementary material to Visser & Fiksen, "Optimal foraging in marine

% ecosystems models: selectivity, profitability and switching"

% Marine Ecology Progress Series 473:XXX?XXX 

% MATLAB code for calculating the grazing rate of a planktonic predator on

% size spectrum of motile and non-motile prey based on 

% 1) prey profitability

% 2) the predators ability to choose between feeding modes

% This code reproduces the switching behaviour illustrated in Figure 6.

% We invite interested readers to explore this behaviour and its

% consequences for different parameter settings.

% The code is less efficient than it could be, but this is done for the

% sake of clarity.

%

Npre = 101;  %number of prey classes to be simulated

Nmot = 21;  % number of motile - nonmotile ratios to be simulated

% initiate matices and vectors

gmotfitclass = zeros(Nmot,Npre);

gnonfitclass = zeros(Nmot,Npre);

gmotmaxclass = zeros(Nmot,Npre);

gnonmaxclass = zeros(Nmot,Npre);

mcruiseinc = zeros(Nmot,Npre);

mambushinc = zeros(Nmot,Npre);

gcruisemax = zeros(Nmot,1);

gambushmax = zeros(Nmot,1);

gambushmax = zeros(Nmot,1);

gfitness = zeros(Nmot,1);

gmaxmode = zeros(Nmot,1);

rcruisemin = zeros(Nmot,1);

rcruisemax = zeros(Nmot,1);

rambushmin = zeros(Nmot,1);

rambushmax = zeros(Nmot,1);

 

% size classes for both motile and nonmotile prey 

i = linspace(-6,-2,Npre);

r = 10.^i       % prey sizes (m)

R = 10^-2;      % predator size (m) 

ccell = 1e5;    % prey cell carbon content gC/m3

vol = 4/3 * pi * r.^3;  %prey cell volume (m^3)

 

% profitability functions

ene = ccell * vol;          % prey cell content (gC) Eq(3)

han = 1 + 1000*(r/R).^3;    % handling time (s) Eq(4)

phi = (1 - r/R).^3;         % catchibility Eq(5)

pro = phi.*ene./han;        % profitability (gC/s) Eq (1)

 

% sort profitability: jsort stores the indices of the size classes in order

% of descending profitability

[psort,jsort] = sort(pro,'descend');    

 

% the relative abundance of motile to nonmotile prey

motile = linspace(0,1,Nmot);

nonmot = 1 - motile;

 

totalbiomass = 1e-1; % total biomass in environment gC/m^3

% biomass in each size class for motile and nonmotile prey. Simplest case 

% assumes a linear spectrum with equal biomass per logarithmic interval.

genrelbiomass = totalbiomass*ones(size(r))/Npre;

    

% Encounter kernel for cruise - ambush feeding mode on motile and

% non-motile prey. Units (m^3/s). Simplest case reflecting symmmetry

% between motile and non-motile targets and searchers. Assumes encounter

% cross section depends on the size of both encounter partners, the (R +

% r)^2 term, and the relative velocity is proportional to the size of the

% moving encounter partner, R or r for cruise and ambush respectively. 

% Encounter rates between stationary patners and co-moving partners set to

% zero but can be readily modified when needed. The speed of swimming is

% assumed a linear function of size, and the overall scaling factor sets

% clearance rates to be consistent with general trends observed by Kiørboe

% 2011.

bnoncruise = 30*R*(R + r).^2;    % 

bnonambush = zeros(size(r));     % 

bmotcruise = zeros(size(r));     % 

bmotambush = 30*r.*(R + r).^2;   %

% As a non-zero motile-cruise encounter example we offer:

% bmotcruise = 5*(1 - r/R).*sqrt(R^2 + r.^2).*(R  + r).^2;     

% where (1 - r/R) indicates escape probability

%       sqrt(R^2 + r.^2) gives relative speed

%       (R  + r).^2 is the encounter cross sectional area as above.

 

% the fitness trade-off

    costcruise = 1e-10; % additional cost for cruise feeding mode gC/s

    mucruise = 4;       % additional mortality risk of cruise feeding. 

% Mortality risk is assumed here to be a constant factor greater than

% mortality risk of being stationary

% Note that cost and mortality are essentially treated as marginal

% increases due to cruise feeding. Thus the marginal cost and risk for

% ambush are by definition zero:

    costambush = 0;     

    muambush = 0;   

    

for i = 1:Nmot % loop over different motile-nonmotile ratios

 

 

% number density per prey class (cells / m^3) for motile and nonmotile

% fractions

    amot = motile(i)*genrelbiomass./ene;

    anon = nonmot(i)*genrelbiomass./ene;

 

% encounter rate for each feeding mode

    abcn = bnoncruise(jsort).*anon(jsort);

    abcm = bmotcruise(jsort).*amot(jsort);

    aban = bnonambush(jsort).*anon(jsort);

    abam = bmotambush(jsort).*amot(jsort);

 

 

% Cruise feeding

% Finding the maximum grazing rate for profitability ranked prey classes is

% relatively straightforward.

    % The numerator of Eq(2) evaluated over sucessive ranges of

    % profitability 

    k1cruise = cumsum((abcn + abcm).*phi(jsort).*ene(jsort)); 

    % Likewise, the contribution to the denominator of Eq(2) evaluated over

    % sucessive ranges of profitability 

    k2cruise = cumsum((abcn + abcm).*han(jsort));               

    % Grazing rate vector including successively less profitable prey Eq(2)

    gcruise = k1cruise./(1 + k2cruise); 

    % Find the maximum grazing rate for cruise feeding mode

    [gcruisemax(i), jcruisemax] = max(gcruise);

    jcruiseinc = psort > gcruisemax(i);

    mcruiseinc(i,:) = pro > gcruisemax(i);

    % Optimal diet breadth

    rcruisemin(i) = min(r(mcruiseinc(i,:)==1));

    rcruisemax(i) = max(r(mcruiseinc(i,:)==1));

    % Grazing pressure per size class

    gcmclass(jsort) = jcruiseinc.*abcm.*phi(jsort).*ene(jsort)...

        /(1 + k2cruise(jcruisemax));

    gcnclass(jsort) = jcruiseinc.*abcn.*phi(jsort).*ene(jsort)...

        /(1 + k2cruise(jcruisemax));

    gcruiseclass = gcmclass + gcnclass;

 

% Ambush feeding : similar to above

    k1ambush = cumsum((aban + abam).*phi(jsort).*ene(jsort));

    k2ambush = cumsum((aban + abam).*han(jsort));

    % Grazing rate vector including successively less profitable prey

    gambush = k1ambush./(1 + k2ambush);

    [gambushmax(i), jambushmax] = max(gambush); % maximum grazing rate

    jambushinc = psort > gambushmax(i);

    mambushinc(i,:) = pro > gambushmax(i);

    % Optimal diet breadth

    rambushmin(i) = min(r(mambushinc(i,:)==1));

    rambushmax(i) = max(r(mambushinc(i,:)==1));

    % Grazing pressure per size class

    gamclass(jsort) = jambushinc.*abam.*phi(jsort).*ene(jsort)...

        /(1 + k2cruise(jambushmax));

    ganclass(jsort) = jambushinc.*aban.*phi(jsort).*ene(jsort)...

        /(1 + k2cruise(jambushmax));

    gambushclass = gamclass + ganclass;

 

% Fitness estimates for feeding modes Eq(13)

    fcruise = (gcruisemax(i) - costcruise)/(1 + mucruise);

    fambush = (gambushmax(i) - costambush)/(1 + muambush);

 

% Deduce optimal (preferred) feeding mode

    if fcruise > fambush

        gmaxmode(i) = gcruisemax(i);

        prefmode = 'cruise';

        gmotmaxclass(i,:) = gcmclass;

        gnonmaxclass(i,:) = gcnclass;

    else

        gmaxmode(i) = gambushmax(i);

        prefmode = 'ambush';

        gmotmaxclass(i,:) = gamclass;

        gnonmaxclass(i,:) = ganclass;

    end

    

% A fitness weighted grazing estimate Eq(14)

    qcruise = fcruise/(fcruise + fambush);

    qambush = fambush/(fcruise + fambush);

    gfitness(i) = qcruise*gcruisemax(i) + qambush*gambushmax(i);

% Grazing pressure per size class for each feeding mode

    gmotfitclass(i,:) = qcruise*gcmclass + qambush*gamclass;

    gnonfitclass(i,:) = qcruise*gcnclass + qambush*ganclass;

 

disp(['motility fraction ' num2str(motile(i))...

      ' preferred feeding mode: ' prefmode ])

 

end

 

% output

figure(1)

clf;

day2sec = 24*3600;

plot(motile,gambushmax*day2sec,'og',motile,gcruisemax*day2sec,'ob'...

    ,motile,gfitness*day2sec,'-r',motile,gmaxmode*day2sec,'-k')

xlabel('motility ratio')

ylabel('grazing rate (gC/day)')

legend('ambush','cruise','fitness','max mode')

 

figure(2)

clf;

semilogy(motile,rcruisemin,'-g',motile,rcruisemax,'--g'...

        ,motile,rambushmin,'-b',motile,rambushmax,'--b')

legend('cruise min','cruise max','ambush min','ambush max')

xlabel('motility ratio')

ylabel('optimal diet breadth')



