Below are exercises and MatLab scripts that support the mathematical biology module of the joint ASU-SCC REU.

**Exercises**

Ebola project

Lemming-fox project

Bear-seal project

Bat rabies project

ITB_PredPrey: This Module comes from a draft of a book I’m preparing called *Introduction to Mathematical Biology*. It shows the phase plane analysis of the model using both modern techniques and the method used by Lotka.

Linear stability analysis, worked example

**MatLab Scripts (programs)**

The following file (LogisticMap) is a MatLab .m file that runs the logistic map from the introductory exercise. It also makes a cobweb diagram. Due to security constraints on this server, it had to be saved as a .doc file instead of MatLab’s standard .m. Therefore, to run the script you’ll need to do the following:

- Download the .doc file but
**do not**open it in Word or any other word processor. - Open the file in a text editor (or even in the MatLab editor).
- Copy the code completely from the text editor.
- Open a new script file in MatLab.
- Paste the code into the new file and save it under any name with a .m extension.

MatLab logistic model script: LogisticMap

The following example scripts can be copied and pasted into the MatLab editor and saved with any file name followed by a .m extension. Beware, though, that you may find formatting problems if you just cut-and-paste, especially the quotation marks.

Experienced programmers may note that these are not the most efficient algorithms, but they have the advantage of transparency, making it easier for someone who is not a MatLab jockey to follow the logic.

% LOGISTIC MODEL TEACHING MODULE

%

% This script runs the discrete-time logistic map.

%

% Written by John Nagy

% June 2013 (revised May 2016)

%———————————————————————-

% Initialize and set parameters

%———————————————————————-

% Clear graphical and work spaces

clf

clear

hold off

% Set demographic parameters (birth, death, competition rates)

lambda = 2;

mu = 0.6; % Must be between 0 and 1, inclusive

c = 0.1;

% Compute the derived parameters

r = lambda – mu;

K = r/c;

% Set the time horizon

maxT = 50;

% Preallocate the main data vector and set the initial number of animals

x = zeros(maxT,1);

x(1) = 5; % In hundreds.

%———————————————————————-

% Run the model

%———————————————————————-

for i=2:maxT

x(i) = x(i-1) + r*x(i-1)*(1-x(i-1)/K);

end

%———————————————————————-

% Plot 1: Time series

%———————————————————————-

plot(1:maxT, x,’-o’) % Plots the basic time series as open points and lines

ax = gca; % Gain control of the current plot axes

line([0 maxT], [K K], ‘color’, ‘red’) % Plots carrying capacity in red

ax.XAxis.FontSize = 12; % Adjust font size of axis marks (numbers)

ax.YAxis.FontSize = 12;

xlabel(‘Time (Number of generations)’, ‘interpreter’, ‘latex’, ‘FontSize’, 14)

ylabel(‘Number of moose (Hundreds)’, ‘interpreter’, ‘latex’, ‘FontSize’, 14)

%———————————————————————-

% Plot 2: Cobweb plot

%———————————————————————-

figure % Create a new plot window

% Calculate the positive portion of the r.h.s of the model (return map)

% [Next 3 lines]

xbaseMax = K*(1+r)/r;

xbase = 0:0.01:xbaseMax;

f = xbase + r*xbase.*(1-xbase/K);

% Plot the return map and the identity line

% [Next 3 lines]

plot(xbase, f, ‘-b’)

ax=gca; % see above

ax.XAxis.FontSize = 12; % see above

ax.YAxis.FontSize = 12;

xlabel(‘Number of moose in generation $n$ (hundreds)’, ‘interpreter’, ‘latex’, ‘FontSize’, 14)

ylabel(‘Number of moose in generation $n+1$ (hundreds)’, ‘interpreter’, ‘latex’, ‘FontSize’, 14)

line([0 ax.XAxis.Limits(2)], [0 ax.XAxis.Limits(2)], ‘color’, ‘black’)

% Overlay the solution x onto the return map

% [Next 2 lines]

hold on

plot(x(1:(maxT-1)),x(2:maxT),’or’)

% Plot the cobweb

line([x(1) x(1)], [0 x(2)], ‘Color’, ‘r’)

line([x(1) x(2)], [x(2) x(2)], ‘Color’, ‘r’)

for i = 2:(maxT-1)

line([x(i) x(i)], [x(i), x(i+1)], ‘Color’, ‘r’)

line([x(i) x(i+1)], [x(i+1) x(i+1)], ‘Color’, ‘r’)

end

The following script requires the function gen_roots, which must be copied verbatim with the filename “gen_roots.m” (including underscore) and saved in the same file as this one. The gen_roots function is listed after this script.

------------------------------------------------------------------

% MCTP example Nullcline Plot % % See Linear Statbility Theory handout % % Written by John Nagy % June 2013

% Parameterize

K = 1; r = 0.2; a = 2;

Max_u = r/a;

% Define the horizontal axis (predators) u = eps:0.001*Max_u:Max_u;

v_null = K*(1 - u*a/r);

u_null = u./(1-exp(-a*u));

plot(u,v_null, u,u_null,'-r')

xlabel('Number of predators (\itu)') ylabel('Number of prey (\itv)')

% Find the actual fixed point % First, subtract the two nullclines:

Null = v_null - u_null;

%Then fine the root:

uFP = gen_roots(u, Null); vFP = K*(1 - uFP*a/r);

% Plot the fixed points

line([uFP uFP], [0 K], 'color', 'black', 'LineStyle', '--') line([0, Max_u], [vFP vFP], 'color', 'black', 'LineStyle', '--')

% In the plot, write the values of the fixed points: % First, convert the values to strings and write the string to be placed % in the plot

uFP = num2str(uFP); vFP = num2str(vFP);

%The next construction concatinates the strings into a single string. uStr = strcat('u fixed point = ', uFP); vStr = strcat('v fixed point = ', vFP);

%Now, write the text on the plot text(0.06, 0.8, uStr) text(0.06, 0.76, vStr)

------------------------------------------------------------------

function rset = gen_roots(d,f) % Returns a column vector, rset, that represents the roots of a general % function, f, obtained by simple linear interpolation It takes two % vectors, d (domain, x-values) and f (function, y-values). If no root % (zero) exists, then rset = NaN. % % Written December 2009 % Written by John D. Nagy

len = length(d); s = sign(f);

if abs(sum(s)) == len rset = NaN; % the abs val of sum of sign function = len iff no roots return end

% Now we know there's at least one root

rset=[];

for i=1:len-1 if f(i) == 0 rset=[rset; d(i)]; elseif s(i)+s(i+1) == 0 interp = (f(i)*d(i+1)-f(i+1)*d(i))/(f(i)-f(i+1)); rset=[rset; interp]; end end