%Geoffrey Johnston
%2-2-09
%Modeling the behavior of mosquito populations under drug pressure

%The distribution is the distance of each of the progeny from the parent.
%We are charting this over time.
%Actually, it is the distribution of expected progeny from the parent.
%We are running different distributions in competition with one another.
%A low variance means that the DNA is stable.
%A high variance means that the DNA is unstable, i.e. is highly mutative.


%user adjustable parameters:
natmutrate = 2*10^-3; %we are looking only for beneficial mutations here...
griddimrow = 5;
griddimcol = 5;
time = 20;
fxdmn = 0;
begmn = 0;
begvar = 1;
immunemn = zeros(griddimrow,griddimcol);
drugpressmin = .1;
drugpressmax = 5;
drugpress = drugpressmin*ones(griddimrow,1);
xmn = -5; %these are the axes for the hist plots at end
xmx = 5;
klngspred = 1; %Damn, the model is sensitive to this parameter!
klngparam = 2;
counter = 1;
cnt = 1;

%Let's assume for now that the places with the most drug pressure have the
%least number of immune individuals... these places can support higher
%parasite loads... that is... drug press high, parasite loads high
%numofind

if griddimrow > 1,
    
for i = 2:griddimrow
    drugpress(i,1) = drugpress(i-1)+ (drugpressmax - drugpressmin)/(griddimrow-1);
end

for i = 2:griddimcol
    drugpress = cat(2,drugpress,drugpress(:,1));
end

end

drugpress2 = drugpressmin*ones(1,griddimcol);

if griddimcol > 1,
    
for i = 2:griddimcol
    drugpress2(1,i) = drugpress2(i-1)+ (drugpressmax - drugpressmin)/(griddimcol-1);
end

for i = 2:griddimrow
    drugpress2 = cat(1,drugpress2,drugpress2(1,:));
end

end

if griddimrow == 1 && griddimcol~=1,
    drugpress = drugpress2;
    
else
    
    if griddimcol == 1 && griddimrow ~=1,
        
        drugpress2 = drugpress;
    else
        
    if griddimcol == 1 && griddimrow ==1,
        
        drugpress = .5*drugpressmax;
        drugpress2 = .5*drugpressmax;
    end
    
    end
    
end

drugpress = .5*drugpress + .5*drugpress2;
numofind = round(1000*drugpress);
maxnumofind = max(max(numofind));

%This ends the user-defined section.


%let's create the right number of species
%follow along with sub2ind(drugpress,griddimrow,griddimcol)


for i = 1:griddimrow*griddimcol;
   
        spec(i,:)= random('norm',begmn,begvar,1,maxnumofind);
    %each row of this matrix contains all of the data... it is horizontal,
    %not vertical
end



for i = 1:griddimrow;
   
    for j = 1:griddimcol;

    subplot(griddimrow, griddimcol, cnt);
    hist(spec(cnt,[1:numofind(i,j)]))
    title('No dp')
    %set(gca,'xtick',[],'ytick',[]);
    xlim([xmn xmx]);
    cnt = cnt+1;

    end
end
cnt = 1;
% we could put this population on a grid and then spead the information,
% i.e. mean and variance, but for now let's just make two pops

% we will need ITERATIONS and pressure; we will let the mean and variance
% vary RANDOMLY!  (note that we could add higher moments later...!)

%LET THE SPECIES NAVIGATE THE SELECTION LANDSCAPE

%now let's apply some drug pressure
% we will kill off the lots of the guys... either the most unfit or the
% tails at one or both ends

%note that our fxdmn, the one targeted by the drug, is 0, ALWAYS! 


for i = 1:griddimrow;
    for j = 1:griddimcol;
    
[specmn(i,j),specvar(i,j)] = normfit(spec(counter,[1:numofind(i,j)]));
[newsp(counter,:),immunemn(i,j),newmn(i,j),newvar(i,j)] = killingfntm2(spec(counter,[1:numofind(i,j)]),time,fxdmn,numofind(i,j),drugpress(i,j),specmn(i,j),specvar(i,j),max(max(drugpress)),klngspred,klngparam,maxnumofind);

counter = counter + 1;

    end

    %still need to pass info to neighboring cells!
end

figure
counter = 1;

for i = 1:griddimrow;
   
    for j = 1:griddimcol;

    subplot(griddimrow, griddimcol, cnt);
    hist(newsp(cnt,[1:numofind(i,j)]))
    xlabel(sprintf('%.3f dp.',drugpress(i,j)));
    xlim([xmn xmx]);
    cnt = cnt+1;

    end
end

cnt = 1;

figure
imagesc((newvar-specvar));
title('This is a plot of the change in variances.');
xlabel('dp increases from right to left.')
ylabel('dp increases from top to bottom.')
colorbar;

figure
if griddimrow>1 && griddimcol>1,
surf((newvar-specvar));
else
imagesc((newvar-specvar));
end

title('This is a plot of the change in variances.');
xlabel('dp increases from right to left.')
ylabel('dp increases from top to bottom.')
colorbar;


figure 
imagesc(drugpress);
title('This is a plot of the drug pressure.')
colorbar;

figure 
imagesc(immunemn);
title('This is a plot of the new immune means.')

figure
imagesc((newmn - specmn));
title('This is a plot of the change in means.')
colorbar;

%now we need to compete and repopulate these guys!
%so they were in the host, now they are outside...
%actually the regrowth should happen in the host...

%the mosquito will then cull the parasites at their fractional levels
%randomly

%genetic drift on the last few...

%Now we have our grid of drug pressure mutants.  These have changed their
%variances in order to complete the most effective linear search.  Let's
%now use these variances to give the mutation rate...
%as the variance increases, the mutation rate increases...

mutrate = zeros(griddimrow,griddimcol);

for i = 1:griddimrow;
    
    for j= 1:griddimcol;
    
        mutrate(i,j) = natmutrate*newvar(i,j);
        
    end
    
end

%Now we will use this mutation rate to predict the number of mutants
%produced... this is the mutation rate times the number of parasites in the
%body... influenced by immune action, relative success rates...

%I assumed that all of the individuals could support the same capacity of
%parasites.... this is not true.  Rather, some people should be able to
%support more than others...

%We will adjust this by changing the numofind into a matrix...
%Now find the number of mutants expected...

nummut = zeros(griddimrow,griddimcol);

for i = 1:griddimrow;
    
    for j= 1:griddimcol;
    
        nummut(i,j) = mutrate(i,j)*numofind(i,j);
        
    end
    
end

%This is the reason that the mutants with higher variance are more fit...
%the mutants with higher var have more mutations on average, and so there
%is more chance that they will survive to the next round!

%Recall that variation in and of itself is not beneficial.
%It is beneficial only insofar as it increases the number of beneficial
%mutations.

%This is where the analytics come in....

%So, in retrospect, why did the mutators win?  Because they could produce
%more mutations...