%%% This is Matlab code for the manuscript
%%% "Snail hosts abundance mediates the effects of antagonist interactions
%%% between trematodes on the transmission of human schistosomes"
%%% by Douchet, Haegeman, Allienne, Boissier, Senghor and Rey.
%%% The code generates Figure 2 of the main text.

% set parameter values
p=struct('np',nan,'g',.2,'e',.01,'fS',3,'K',nan, ...
    's',2,'u',1,'dM',.01,'dI',.1,'dC',.1);

% compute equilibrium values
Ns=1:6; % values for number of trematode species
Ks=[1 10 100]; % values for carrying capacity
for cK=1:length(Ks)
    p.K=Ks(cK);
    res=zeros(length(Ns),3);
    for cN=1:length(Ns)
        p.np=Ns(cN);
        if cN==1
            x0=.1*ones(1,4);
        else
            x0=.1*ones(1,p.np+4);
        end
        opts=odeset('nonnegative',1:length(x0));
        % simulate during long time, here 1000 time units
        [T,X]=ode45(@(t,x) divpara_rhs(x,p),[0 1000],x0,opts);
        % extract equilibrium values
        if cN==1
            res(cN,:)=[0 X(end,3)/sum(X(end,2:3)) 1];
        else
            res(cN,:)=[sum(X(end,4:end-2))/sum(X(end,3:end-2)) ...
                sum(X(end,3:end-2).*(1:cN)/cN)/sum(X(end,2:end-2)) ...
                X(end,end-1)/sum(X(end,end-1:end))];
        end
    end
    ress{cK}=res;
end

% plot simulation results
cols=get(gca,'colororder');
for cK=1:length(Ks)
    subplot(1,3,cK)
    plot(Ns,ress{cK}(:,1),'.-', ...
        'color',cols(3,:),'markersize',12,'linewidth',1)
    hold on
    plot(Ns,ress{cK}(:,2),'.-', ...
        'color',cols(2,:),'markersize',12,'linewidth',1)
    plot(Ns,ress{cK}(:,3),'.-', ...
        'color',cols(4,:),'markersize',12,'linewidth',1)
    plot(Ns,1./Ns,'--','color',cols(4,:))
    xlim([1 6])
    xlabel('Number of trematode species')
    ylabel('Equilibrium ratio')
    title(sprintf('Carrying capacity K=%d',Ks(cK)),'fontweight','normal')
end
hl=legend('co-infestation ratio', ...
    'Schisto prevalence','Schisto cercariae ratio');
set(hl,'box','off','location','northeast')
