function sense = periodicSense(P,t,h,n0,cycle,Opt)
% This function calculates the transient sensitivity of the total
% population size at time t after an initial perturbation in the elements
% of the phase h (i.e. P{cycle(h)})
%
% P = an array of cells, where each cell contains a transition matrix
% for one phase of a cycle (the matrices in P can be different
% sizes provided they tally correctly for the specified cycle)
% t = the time after the initial perturbation
% h = the sensitivities for a perturbation in phase h
% n0 = the initial population vector
% cycle = is a vector giving the cyclical order of each phase.
% The default is [1,2,3, ...,H], where H=number of phases
% e.g. cycle = [1 3 2 4] means that the order of the
% transition matrices in a cycle is P{1} followed by P{3}
% and then P{2} and P{4}, so that the transition matrix for the
% whole cycle is P{4}*P{2}*P{3}*P{1}
% Opt = a structure that gives different output options
% Opt.elasticity = true Calculate elasticities
% Opt.immediate = true Calculate immediate sensitivities
% Opt.timescaled = true Calculate time-scaled sensitivities
% [default for all is false]
%
% OUTPUT
% sense(i,j) = the transient sensitivity of the total population size a
% time t after a perturbation of element P{cycle(h)}(i,j)
% Depending on output options this could be within-cycle,
% elasticity, etc
% Written by Jon Yearsley and Shana Mertens, Jan 2006
% jonathan.yearsley@unil.ch
% shana.mertens@bbsrc.ac.uk
if nargin<4,
% Check that there are the right number of inputs
error('periodicSense requires at least 4 input arguments: P, t, h and n0')
elseif nargin<5,
nPhases = length(P); % Number of phases
cycle = 1:nPhases; % Set default cycle
Opt.elasticity = false; % Set default options
Opt.immediate = false;
Opt.timescaled = false;
elseif nargin<6,
nPhases = length(cycle);
Opt.elasticity = false; % Set default options
Opt.immediate = false;
Opt.timescaled = false;
else
nPhases = length(cycle);
% Check the options are set properly
if ~isfield(Opt,'elasticity')
warning('Setting elasticity option to default')
Opt.elasticity = false;
end
if ~isfield(Opt,'immediate')
warning('Setting immediate option to default')
Opt.immediate = false;
end
if ~isfield(Opt,'timescaled')
warning('Setting timescaled option to default')
Opt.timescaled = false;
end
end
if h>length(cycle),
% Check that the required phaseis consistent with the cycle information
error(['The argument h is too large. There are only ' ...
num2str(nPhases) ' in the current cycle'])
end
nStages = zeros(length(P),2);
for p=1:length(P),
nStages(p,:) = size(P{p}); % Number of stages in each transition matrix
end
% Finally check that the cycle is possible given the dimension of the
% matrices in P
if any(nStages(cycle,1) ~= nStages(cycle([2:nPhases,1]),2)),
error('The matrix dimensions of P are wrong for the given cyclic permutation')
end
% Initialise the matrix to contain the sensitivities
sense = zeros(nStages(h,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate decomposition of sensitivity
%
% Now calculate D^{(+h)} and D^{(-h)} (equations 9 and 10 in the paper)
% which when multiplied together give Dh
% Dh multiplies the "cycle" sensitivities and gives us
% the sensitivity of a particular element in phase h
Dminus_h = eye(nStages(cycle(end),1));
for i=1:h-1,
Dminus_h = P{cycle(i)} * Dminus_h;
end
Dplus_h = eye(nStages(cycle(h),1));
for i=h+1:nPhases,
Dplus_h = P{cycle(i)} * Dplus_h;
end
% Finally create the transition matrix for an entire cycle
A = Dplus_h * P{cycle(h)} * Dminus_h;
% Calculate the deviation of n0 from the stable stage structure
% Calculate right eigenvectors and eigenvalues
optEig.disp = 0;
[w1,~] = eigs(A,1,'lm',optEig);
w1= abs(w1); % Make sure sign of w1 is positive
% The equilibrium stage distribution
if length(n0)~=length(A),
error('The number of stages in n0 does not agree with the model')
end
n0 = reshape(n0,length(A),1); % This ensures n0 is a column vector
n_infty = sum(n0)/sum(w1); % Use unstandardised right eigenvector because
% S_h is calculated with unstandardised
% eigenvectors
% Deviation of initial population structure from asymptotic stable stage
% structure
delta_n = n0 - sum(n0)/sum(w1)*w1;
if Opt.elasticity,
%Calculating population size at time t
n = A^t * n0;
nTot=sum(n); % Total population size
if Opt.immediate,
n_within = Dplus_h \ n;
nTot = sum(n_within);
end
end
for k=1:nStages(h,1),
for l=1:nStages(h,2),
% This is equation 8 in the Ecology paper
Dh = Dplus_h(:,k) * Dminus_h(l,:);
% Find the elements of Dh which are non-zero
[iInd,jInd] = find(Dh);
% Now do the same thing only break up the sensitivity into two parts, a la
% Ecological modelling paper (Yearsley, 2004)
S_h = zeros(length(n0),1);
deltaS_h = zeros(length(n0));
if ~isempty(iInd),
for i=1:length(iInd);
S_h = S_h + Dh(iInd(i),jInd(i)) * real(S(A,iInd(i),jInd(i),t));
deltaS_h = deltaS_h + Dh(iInd(i),jInd(i)) * real(delta_S(A,iInd(i),jInd(i),t));
end
% And finally calculate the sensitivity of the population size
% based upon the decomposition in Ecological Modelling paper
% (there are other ways to do this without using a
% decomposition)
% Note that we're claculating the sensitivity of total
% population size, hence the sum( * ,1)
if Opt.immediate,
% sense(k,l) = sum(Dplus_inv * (S_h * n_infty + deltaS_h * delta_n),1);
sense(k,l) = sum(Dplus_h \ (S_h * n_infty + deltaS_h * delta_n),1);
else
sense(k,l) = sum(S_h * n_infty + deltaS_h * delta_n,1);
end
end
end
end
if Opt.elasticity,
% And elasticities, scaled for time
sense = P{cycle(h)} .* sense / nTot;
end
if Opt.timescaled
sense = sense / t;
end