Alright, thanks Roger! I've included the entirety of the code here.
%(*this is the code for the branching tree algorithm for cancer mutations*)
%(*this is where we define our parameters*)
b = 0.25; %(*birth rate*)
s = 10^2; %(*decrease in death rate for more fit cells*)
uA = 10^3; %(*mutation rate for advantageous mutants*)
uN = 10^2; %(*mutation rate for neutral mutants*)
topTime = 600; %(*number of time steps for each run*)
posMut = 30; %(*length of our mutation vector*)
NumRuns = 10000; %(*number of runs*)
deaths = zeros(posMut);
bplusdeaths = zeros(posMut);
for j=1:posMut,
deaths(1,j) = b*(1  s);
end
for i=2:posMut,
for j=1:posMut,
deaths(i,j)=deaths(i  1,1)*(1  s);
end
end
for j=1:posMut,
for i=1:posMut,
bplusdeaths(i,j) = b + deaths(i,j);
end
end
totalAbundance = zeros(posMut);
averageAbundance = zeros(posMut);
goodRunsCounter = 0;
%(*here is where we loop through each of the runs*)
for q = 1:NumRuns,
abundance = zeros(posMut);
abundance(1,1) = 1;
time = 0;
%(*find gamma*)
while(time<topTime)
gamma = 0;
for i=1:posMut,
gamma = gamma + dot(abundance(i,:),bplusdeaths(i,:));
end
if(gamma == 0)
break;
end
%(*Exponential distribution*)
deltaT = log(rand)/(gamma);
handOfGodCol = rand(1);
relAbundCol = cumsum(abundance);
for i=2:posMut,
relAbundCol(posMut,i) = relAbundCol(posMut,i) + relAbundCol(posMut,i1);
end
relAbundOne = relAbundCol(posMut,:)/sum(sum(abundance));
j=1;
chosenMutantN = 1;
%(*this chooses which neutral mutant we are choosing to have something happen to*)
while(relAbundOne(j)<=handOfGodCol)
chosenMutantN = j + 1;
j = j + 1;
end
handOfGodRow = rand(1);
relAbundRow = zeros(posMut);
relAbundRow(1,chosenMutantN) = abundance(1,chosenMutantN);
for i=2:posMut,
relAbundRow(i,chosenMutantN) = abundance(i,chosenMutantN) + relAbundRow(i1,chosenMutantN);
end
relAbundTwo = relAbundRow(:,chosenMutantN)/sum(abundance(:,chosenMutantN));
j=1;
chosenMutantA = 1;
%(*this chooses which advantageous mutant we are choosing to have something happen*)
while(relAbundTwo(j)<=handOfGodRow)
chosenMutantA = j + 1;
j = j + 1;
end
handOfGodAgain = rand(1);
%(*this chooses what happens to that chosen mutant*)
total = abundance(chosenMutantA,chosenMutantN)*bplusdeaths(chosenMutantA,1);
birth = abundance(chosenMutantA,chosenMutantN)*b*(1  uA  uN)/total;
mutationN = abundance(chosenMutantA,chosenMutantN)*b*(1  uA)/total;
mutationA = abundance(chosenMutantA,chosenMutantN)*b/total;
if(handOfGodAgain < birth)
abundance(chosenMutantA,chosenMutantN) = abundance(chosenMutantA,chosenMutantN) + 1;
elseif(handOfGodAgain < mutationN)
abundance(chosenMutantA,chosenMutantN + 1) = abundance(chosenMutantA,chosenMutantN + 1) + 1;
elseif(handOfGodAgain < mutationA)
abundance(chosenMutantA + 1,chosenMutantN) = abundance(chosenMutantA + 1,chosenMutantN) + 1;
else
abundance(chosenMutantA,chosenMutantN) = abundance(chosenMutantA,chosenMutantN)  1;
end
time = time + deltaT;
end
%(*this is where we find the average outputs of the runs*)
if(gamma ~= 0)
goodRunsCounter = goodRunsCounter + 1;
disp('we just got a good run!');
disp(q);
for v = 1:posMut,
for w=1:posMut,
totalAbundance(v,w) = totalAbundance(v,w) + abundance(v,w);
end
end
end
end
disp('Total Abundance Matrix = ');
disp(totalAbundance);
disp('Number of Runs with surviving cells = ');
disp(goodRunsCounter);
for g=1:posMut,
for h=1:posMut,
if(goodRunsCounter~=0)
averageAbundance(g,h)=totalAbundance(g,h)/goodRunsCounter;
end
end
end
disp('Average Abundance Matrix = ');
disp(averageAbundance);
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gjmogm$evm$1@fred.mathworks.com>...
> "Nick StuckyMack" <njstucky@fas.harvard.edu> wrote in message <gjmk9m$5ke$1@fred.mathworks.com>...
> > Below is an excerpt of my code. I keep getting the error message "??? Attempted to access relAbundCol(NaN,:); index must be a positive integer or logical." The weird thing is that this only shows up after it has run successfully multiple times (the code below is part of a larger loop). I can't seem to figure out what is causing it to come up as not a number. Help!
> > Relevant excerpt of the code is below.
> >
> > b = 0.25; %(*birth rate*)
> > s = 10^2; %(*decrease in death rate for more fit cells*)
> > uA = 10^3; %(*mutation rate for advantageous mutants*)
> > uN = 10^2; %(*mutation rate for neutral mutants*)
> > topTime = 600; %(*number of time steps for each run*)
> > posMut = 30; %(*length of our mutation vector*)
> > NumRuns = 10000; %(*number of runs*)
> >
> > for q = 1:NumRuns,
> > abundance = zeros(posMut);
> > abundance(1,1) = 1;
> > time = 0;
> >
> > relAbundCol = cumsum(abundance);
> >
> > for i=2:posMut,
> > relAbundCol(posMut,i) = relAbundCol(posMut,i) + relAbundCol(posMut,i1);
> > end
> >
> > relAbundOne = relAbundCol(posMut,:)/sum(sum(abundance));
>
> I would advise you to show more of the significant parts of your code if you want to receive useful answers rather than wild guesses. There is very little that can be deduced from what you show here.
>
> Roger Stafford
>
