% Modeling the phenomenon of spruce budworm epidemics
% With a focus on the manipulation of critical points as a means of
% controlling the epidemic.

% By Justin Dunmyre
% Modified by Gavin LaRose, 9 Aug 2012

% Description
% When spruce budworms (not actually a type of worm!) reach epidemic
% levels, they have been stated to be the 2nd largest threat to forests,
% behind forest fires.
%
% This demo gives a reduced model of budworm population (see subfunction
% budwormODE) and with r=.5 and q=10, we see sustained elevated levels of
% the budworm density, u.  
% 
% Accepting the interpretation of r as a birth rate (adjusted for predation
% by birds), and q as a combination of available foliage for consumption and
% threshold for predation, we can ask the question:
%     What can we do to prevent an outbreak of spruce budworms?
% A reasonable first solution is to use insecticide, which would correspond to
% lowering the value of q (because we are reducing the available foliage for
% consumption).   
%
% Pursuing this, the demo first plots the solution to the ODE for various
% levels of q (first individually, then all together for comparison.
% We find that for q at and below 7.3878, starting with elevated population 
% u(0)=3, the population drops to a "refuge" level of approximately u=0.8.  
%
% To explain this phenomenon, we graph the right hand side of the ODE to
% show the critical points (RHS=0).  We see that for q below 7.3878, two of the 
% critical points disappear, leaving a critical point at approximately u=0.8
% This critical point is stable (i.e. solutions will approach it), and 
% corresponds to a density of budworms below epidemic levels.  
%
% Finally, we show the relationship between critical points and the solution
% curve, by simultaneously animating the graph of the solution vs. time, and
% the graph of the solution on the phase line.  
% 
% Insecticide is not always desirable, especially in areas designated for
% hiking and camping.  There are  budworm specific diseases, the
% introduction of which could change the value of r.  Could the introduction
% of such a disease prevent an outbreak?

% (c)2012 Regents of the University of Michigan
% This work is licensed under the
% Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% or send a letter to
% Creative Commons,
% 444 Castro Street, Suite 900,
% Mountain View, CA, 94041, USA.

function budworm()

% Stuff that is easy for you to change
    % Various values of q to use, see budwormODE below
    interestingQ=[10 9 8 7.3878 7 6];
    % These are the values of q for which we will animate solutions
    animatedIndices=[1,3,5];

    % Initial condition for u
    ic = 3;

    % colors to use for various q values
    colorvec={'w','m','y','c','g','b','k'};

    % r parameter, see budwormODE below.  If you change this, you
    %    probably need to use different q
    r=0.5;
   
    function uprime=budwormODE(t,u)
        % A model representing the potential for spruce budworm epidemic
        % in North American forests, taken from the text "Mathematical
        % Biology" by J.D.Murray, third edition.  The parameters in the 
        % model can be adjusted to fit population data very closely.  
        %---------------------------------------------
        % Notice that the equation is a combination of a logistic equation
        % r u(1-u/q) and the subtraction of a term that is always positive
        % that represents predation by an avian species.
        %
        % u represents the scaled budworm density.  For the parameter values
        % chosen, high values of u (say, u>4) correspond to epidemic levels of
        % budworms. 
        %
        % r represents the combination of the effects of budworm birth rate and
        % predation by avian species
        %
        % q is a carrying capacity, limiting the growth of budworms on the basis
        % of food (eg: total surface area of branches) as well as the threshold 
        % for which predation is well supported

        % Values of the parameters corresponding to actual outbreak data 
        % can be found in the original paper by Ludwig,Jones, and Holling,
        % "Qualitative Analysis of Insect Outbreak Systems:
        % The Spruce Budworm and Forest" published in 1978.
        % The parameters have been (drastically) changed here to highlight the 
        % equilibria for discussion
        %--------------------------------------------

        % For q>7.3878, there is a stable equilibrium at epidemic levels.

        uprime=r.*u.*(1-u/q)-((u.^2)./(1+u.^2));
    end
    
    %----------------------------------------------------------------------
    % Below here changes to the code should not be necessary, unless
    % you want to explore changing r instead of q
    % To make that change, the simplest way is to do a find/replace to
    % change q=interestingQ to r=interestingR, and create the interestingR
    % vector similar to the interestingQ vector above.

    % It will be necessary to change the legends and annotations to reflect
    % this change. 
    
    %-----------------------------------
    % First, plot the solution graphs for different values of Q

    scrsz = get(0,'ScreenSize'); % Get the screen size
    pause on
    figh=figure('Position',[1 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
    axes('Parent', figh, 'Color',[0.5 0.5 0.5]);

    hold on;
    for j=1:length(interestingQ)
        q=interestingQ(j);
        [t,u]=ode45(@budwormODE,[0 200], ic);
        line=plot(t,u,[colorvec{j} '-'],'LineWidth',2);
	axis( [ 0 200 0 max(interestingQ) ] );
        ylab = ylabel('u','FontSize',16);
        xlab = xlabel('t','FontSize',16);
        text=annotation(figh,'textbox', [0.4 0.75 0.55 0.15],...
            'Interpreter','none',...
            'String',{'Solution to u''= r u(1-u/q) - (u^2)/(1+u^2)',...
		      ['r=' num2str(r), '; q=', num2str(q)]},...
            'FontSize',20,...
            'FitBoxToText','on');
	pause;
	delete(text);
	delete(line);
	delete(xlab);
	delete(ylab);
    end
    % close(figh);
    hold off;
    clf;

    %-----------------------------------
    % Next, plot all of the solution graphs together

    pause on
    % figh=figure('Position',[1 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
    axes('Parent', figh, 'Color',[0.5 0.5 0.5]);
    hold on;
    ylabel('u','FontSize',16);
    xlabel('t','FontSize',16);
    for j=1:length(interestingQ)
        q=interestingQ(j);
        legendCell{j}=['q = ' num2str(q)]; % Building up the legend
        
        [t,u]=ode45(@budwormODE,[0 200], ic); % solve the ODE
        
        plot(t,u,[colorvec{j} '-'],'LineWidth',2) % plot with the right color
    end
    text=annotation(figh,'textbox', [0.36 0.33 0.53 0.20],...
        'Interpreter','none',...
        'String',{'Solution to u''= r u(1-u/q) - (u^2)/(1+u^2)',
                  ['r=' num2str(r)]},...
        'FontSize',20,...
        'FitBoxToText','on');
    legend(legendCell);
    
    pause
    
    %-----------------------------------
    % Now go and plot the RHS of the ODE against u

    % we'll use this x to discover equilibria
    x=linspace(0,max(interestingQ)+1,10000);
   
    % The next few commands set up the figure window
    figh=figure('Position',[1 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
    axes('Parent', figh, 'Color',[0.5, 0.5, 0.5]);
    hold on; % We'll be plotting multiple curves in the same figure   
    for j=1:length(interestingQ)
        q=interestingQ(j);
        legendCell{j}=['q = ' num2str(q)]; %Building up the legend
        bprime=budwormODE(1,x);

        plot(x,bprime,[colorvec{j} '-'],'LineWidth',2)
    end
    
    % also plot a dotted line on the x-axis
    plot(x,zeros(1,length(x)),'k:','LineWidth',2);

    % add the legend and annotation
    legend(legendCell,'FontSize',16);
    text=annotation(figh, 'textbox', [0.20 0.80 0.51 0.12],...
        'Interpreter','none',...
        'String', {['graph of r u(1-u/q) - (u^2)/(1+u^2)'],...
        ['for various values of q']},...
        'FontSize',20,...
        'FitBoxToText','on');

    axis([0 max(interestingQ)+1 -.5 .5]) %set the axis limits
    ylabel('u\prime');
    xlabel('u');
    pause   

    %-----------------------------------
    % Finally generate animations showing the generation of the
    %    solution curve and corresponding values on the phase line
    
    for k=1:length(animatedIndices)
        figh=figure('Position',[scrsz(3)/4 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
        subplot1=subplot(2,1,1, 'Parent', figh, 'Color',[0.5, 0.5, 0.5]);
        subplot2=subplot(2,1,2, 'Parent', figh, 'Color',[0.5, 0.5, 0.5]);
        hold(subplot1);
        hold(subplot2);

        q=interestingQ(animatedIndices(k));
        [t,u]=ode45(@budwormODE,[0 200], ic);

	% the solution curve is graphed in the top figure
	plot( [t(1), t(1)+.01], [u(1), u(1)+.01*((u(2)-u(1))/(t(2)-t(1)))],...
            [colorvec{animatedIndices(k)} '-'],...
	    'LineWidth', 2, 'Parent',subplot1);
        legend(subplot1, {'Solution to ODE'}, 'location','SE')

        p1 = plot( t(1), u(1), [colorvec{animatedIndices(k)} '.'],...
	    'MarkerSize', 30, 'Parent',subplot1);
        plot(t,u,[colorvec{animatedIndices(k)} ':'],...
            'LineWidth', 2, 'Parent', subplot1) 
        xlabel('t', 'Parent', subplot1);
        ylabel('u', 'Parent', subplot1);

	% then we graph the phase plane in the lower figure
        plot(x,budwormODE(1,x),'k-','LineWidth',2,'Parent',subplot2);
	plot( [u(1), u(1)+.1], [0, 0], [colorvec{animatedIndices(k)} '-'],...
	      'LineWidth', 2, 'Parent',subplot2);
        legend(subplot2,{'RHS of ODE','Solution to ODE'})

        p2 = plot( u(1), 0, [colorvec{animatedIndices(k)} '.'],...
	      'MarkerSize', 30, 'Parent',subplot2);

        % put in the axis as a dotted line
        plot(x,zeros(1,length(x)),'k:','LineWidth',2,'Parent',subplot2);
        
        xlabel('u');
            
        legend(subplot2,{'RHS of ODE','Solution to ODE'})
        axis([0 max(x) -.375 .375]);
        annotation(figh,'textbox', [0.54 0.70 0.17 0.09],...
            'String',{['q=' num2str(q)]},'FontSize',20,...
            'FitBoxToText','on');

        PhaseLineArrowsu=[.25, 1.3, 3.5, 9];
        for j=1:length(PhaseLineArrowsu)
            myu=PhaseLineArrowsu(j);
            if(budwormODE(1,myu)>0)
                PhaseArrow(j)=plot(myu,0,'k>','MarkerSize',15,'MarkerFaceColor','k');
            else
                PhaseArrow(j)=plot(myu,0,'k<','MarkerSize',15,'MarkerFaceColor','k');
            end
        end

        % add the equilibria
        eqn = [ num2str(r), '*x*(1 - x/', num2str(q), ...
            ') - ((x^2)/(1 + x^2))' ];
        evals = solve( eqn );
        for j=1:length(evals)
	    plot( real(evals(j)), 0, 'k.', 'MarkerSize', 30, ...
	        'MarkerFaceColor', 'k', 'Parent', subplot2 );
        end

        pause

	% delete the points showing the initial conditions
	delete(p1);
	delete(p2);

        % do the animation by plotting an increasing percentage of
        %   the solution vector appropriately
        for j=1:length(t)
            tvec=t(1:j);
            uvec=u(1:j);
            solutionDotX=t(j);
            solutionDotY=u(j);

            phaseVecX=u(1:j);
            phaseVecY=zeros(1,j);
            phaseDotX=u(j);
            phaseDotY=0;

            plot( tvec, uvec, [colorvec{animatedIndices(k)} '-'],...
                'LineWidth', 2, 'Parent', subplot1 );
            soldot = plot( solutionDotX, solutionDotY,...
                [colorvec{animatedIndices(k)} '.'], 'MarkerSize',30,...
                'Parent',subplot1 );

            plot( phaseVecX, phaseVecY, [colorvec{animatedIndices(k)} '-'],...
                'LineWidth', 2, 'Parent', subplot2);
            phasedot = plot( phaseDotX, phaseDotY, ...
                [colorvec{animatedIndices(k)} '.'], 'MarkerSize', 30,...
                'Parent', subplot2 );
            
            getframe();
            delete(soldot);
            delete(phasedot);
        end
        pause
        close(figh)
    end
end

% all done
%-------------------------------------------------------------------------------
