########################################## ## model without allelopathic chemicals ## ########################################## library(deSolve) bact.dynamics<-function(t,y,params) { with(as.list(params), { ## initial resource concentrations R1<-y[1] R2<-y[2] ## vector of density of each phenotype of the two species N1<-y[(length(R1)+length(R2)+1):(length(R1)+length(R2)+length(N1))] N2<-y[(length(R1)+length(R2)+length(N1)+1):(length(R1)+length(R2)+length(N1)+length(N2))] #### if N falls below 0.01, call it 0 tmp<-which((N1)<0.01,arr.ind=T) tmp tmp2<-which((N1)>=0.01,arr.ind=T) tmp2 N1b<-array(NA, dim = length(N1)) N1b[tmp]<-(0) N1b[tmp2]<-(N1[tmp2]) N1<-N1b tmp<-which((N2)<0.01,arr.ind=T) tmp2<-which((N2)>=0.01,arr.ind=T) N2b<-array(NA, dim = length(N2)) N2b[tmp]<-(0) N2b[tmp2]<-(N2[tmp2]) N2<-N2b #### change in resource 1 is the sum of resources in through influx and removed by efflux minus those used by sp 1 and sp 2 #### vmax is a funtion of the trait X which has distribution traded-off with f (affinity for resource) and g (nATP produced in metabolism) from Gudelj et al 2008 #### sp1 and sp2 have different shaped tradeoff functions which are defined later dR1.dt<- dil.rate*(R1.food-R1) - sum((1/growth.yield.1)*N1*X*R1/((f1*X)+R1)) - sum((1/growth.yield.1)*N2*X*R1/((f2*X)+R1)) #### R2 is potentially produced and used by both species #### R2 is not brought in through influx but is lost in efflux #### growth on R2 traded-off against growth on R1 dR2.dt<- dil.rate*(-R2)+ sum((1/growth.yield.2)*factor.R2*N1*X*R1/((f1*X)+R1)) + sum((1/growth.yield.2)*factor.R2*N2*X*R1/((f2*X)+R1)) - sum((1/growth.yield.1)*N2*(1-X)*R2/((f2*(1-X))+R2)) - sum((1/growth.yield.1)*N1*(1-X)*R2/((f1*(1-X))+R2)) #### the growth of each spp will be limited by resource abundance so use resources/number of phenotypes in the later equations R1x<-c(R1/length(N1)) R2x<-c(R2/length(N2)) #### The abundance of each phenotype is found by calculating growth on R1 plus that on R2 minus mortality caused by efflux, allelopatic chemicals and stress #### plus the number of individuals that are added to that phenotype by mutation of others (see Gudelj et al 2008) #### The susceptibility of a sp to allelopathic chemicals produced by others is traded off against the tolerance of that sp to stress dN1.dt<-array(NA, dim = length(N1)) dN2.dt<-array(NA, dim = length(N1)) # evolution # growth of N1 using R1 # growth of N1 using R2 # sp1 removed by efflux and stress (S1 is tolerance to stress) dN1a.dt<- mut.rate*(N1[2]-N1[1])+ c.conv*N1[1]*g1[1]*X[1]*(X[1]*R1x/((f1[1]*X[1])+R1x)) + c.conv*N1[1]*g1[1]*(1-X[1])*((1-X[1])*R2x/((f1[1]*(1-X[1]))+R2x)) - dil.rate*N1[1]/S1 dN1.dt[1]<-dN1a.dt for ( i in (2:(length(N1)-1))){ dN1i.dt<- mut.rate*(0.5*N1[i-1]+0.5*N1[i+1]-N1[i])+c.conv*N1[i]*g1[i]*X[i]*(X[i]*R1x/((f1[i]*X[i])+R1x)) + c.conv*N1[i]*g1[i]*(1-X[i])*((1-X[i])*R2x/((f1[i]*(1-X[i]))+R2x) ) -dil.rate*N1[i]/S1 dN1.dt[i]<-dN1i.dt } dN1n.dt<- mut.rate*(N1[length(N1)-1]-N1[length(N1)]) +c.conv*N1[length(N1)]*g1[length(N1)]*X[length(N1)]*(X[length(N1)]*R1x/((f1[length(N1)]*X[length(N1)])+R1x)) + c.conv*N1[length(N1)]*g1[length(N1)]*(1-X[length(N1)])*((1-X[length(N1)])*R2x/((f1[length(N1)]*(1-X[length(N1)]))+R2x)) - dil.rate*N1[length(N1)]/S1 dN1.dt[length(N1)]<-dN1n.dt #### N2 dN2a.dt<- mut.rate*(N2[2]-N2[1])+ c.conv*N2[1]*g2[1]*(1-X[1])*((1-X[1])*R2x/((f2[1]*(1-X[1]))+R2x)) +c.conv*N2[1]*g2[1]*X[1]*(X[1]*R1x/((f2[1]*X[1])+R1x)) - dil.rate*N2[1]/S2 dN2.dt[1]<-dN2a.dt for ( i in (2:(length(N2)-1))){ dN2i.dt<- mut.rate*(0.5*N2[i-1]+0.5*N2[i+1]-N2[i])+c.conv*N2[i]*g2[i]*(1-X[i])*((1-X[i])*R2x/((f2[i]*(1-X[i]))+R2x)) +c.conv*N2[i]*g2[i]*X[i]*(X[i]*R1x/((f2[i]*X[i])+R1x)) -dil.rate*N2[i]/S2 dN2.dt[i]<-dN2i.dt } dN2n.dt<- mut.rate*(N2[length(N2)-1]-N2[length(N2)]) +c.conv*N2[length(N2)]*g2[length(N2)]*(1-X[length(N2)])*((1-X[length(N2)])*R2x/((f2[length(N2)]*(1-X[length(N2)]))+R2x)) +c.conv*N2[length(N2)]*g2[length(N2)]*X[length(N2)]*(X[length(N2)]*R1x/((f2[length(N2)]*X[length(N2)])+R1x)) - dil.rate*N2[length(N2)]/S2 dN2.dt[length(N2)]<-dN2n.dt #### when R falls below 0.01, call it 0 tmp<-which((R1)<0.01,arr.ind=T) tmp2<-which((R1)>=0.01,arr.ind=T) R1b<-array(NA, dim = length(R1)) R1b[tmp]<-(0) R1b[tmp2]<-(R1[tmp2]) R1<-R1b tmp<-which((R2)<0.01,arr.ind=T) tmp2<-which((R2)>=0.01,arr.ind=T) R2b<-array(NA, dim = length(R2)) R2b[tmp]<-(0) R2b[tmp2]<-(R2[tmp2]) R2<-R2b return (list(c(dR1.dt, dR2.dt,dN1.dt,dN2.dt), NULL)) }) } #### specify parameters #dil.rate=0.05 ## rate of addition of resources and removal of resources, individuals and chemicals from the chemostat #eff=0.01 ## number of As produced by each N #R1.food=1000 ## concentration of resources in the influx #c.conv=0.75 ## conversion constant #factor.R2=1 ## how many R2s are produced per R1 used #mut.rate=0.05 ## mutation rate #f1=f1 ## vector describing the affinity of each phenotype of sp 1 for the resource #f2=f2 ## vector describing the affinity of each phenotype of sp 2 for the resource #g1=g1 ## vector describing the number of ATPs produced by each phenotype of sp 1 when using resource #g2=g2 ## vector describing the number of ATPs produced by each phenotype of sp 2 when using resource #X=X ## trait distribution #growth.yield.1=1 ## number of Ns produced by using R1 #growth.yield.2=1 ## number of Ns produced by using R2 #S1=S1 ## tolerance of N1 to stress (defined later). The inverse of S1 is the susceptibility of N1 to allelopathic chemical A2. #S2=S2 ## each phenotype has a given value of X which affects their vmax on R1 (vmax on R2 is a function of 1-X) ## species evolve by change in abundance of phenotypes X<-c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9) f1<-rep(0.1,9) f2<-rep(0.1,9) ## monod shaped tradeoff curves for g g1<-(X/(0.25+X)) plot(g1) g2<-1-(X/(0.75+X)) plot(g2) plot((g1)~X,col=1,pch=16,ylim=c(0,2),type="l") lines(((g2))~X,pch=16,col=2) ## loop from here # data frame for results results<-matrix(NA,nrow=9,ncol=26) for (i in (1:9)){ # stress input values S1<-(i*0.1) S2<-0.9 params <- c(dil.rate=0.05,eff=0.01 ,R1.food=1000, c.conv=0.75, factor.R2=1, mut.rate=0.05,f1=f1,f2=f2,g1=g1,g2=g2,X=X,growth.yield.1=1,growth.yield.2=1,S1=S1,S2=S2) ## make initial populations of N1 and N2 v<-rnorm(9,2.5,sd=1) v<-ifelse(v<0,0,v) V<-sort(v) N1<-c(V[1],V[3],V[5],V[7],V[9],V[8],V[6],V[4],V[2]) v<-rnorm(9,2.5,sd=1) v<-ifelse(v<0,0,v) V<-sort(v) N2<-c(V[1],V[3],V[5],V[7],V[9],V[8],V[6],V[4],V[2]) # Initial values of resources and bacteria y <- c(r1=500,r2=0,n1=N1,n2=N2) # iterations times <- seq(1,1000,0.1) # Solving r <- lsoda(y, times, bact.dynamics, params,atol=0.0001,rtol=0.0001) # put results from final iteration into dataframe. Order: i*0.1, iteration number, abundance of R1, abundance of R2, abundance of sp1 when X=0.1, when X=0.2, when X=0.3 etc, abundance of sp2 when X=0.1, X=0.2 etc results[i,1:22]<-c(i*0.1,as.vector(r[length(r[,1]),])) namefile<-c('c:/folder1/file1.pdf','c:/folder1/file2.pdf','c:/folder1/file3.pdf','c:/folder1/file4.pdf','c:/folder1/file5.pdf','c:/folder1/file6.pdf','c:/folder1/file7.pdf','c:/folder1/file8.pdf','c:/folder1/file9.pdf') pdf(file=namefile[i], height =6, width = 6) matplot(r[,1], r[,2:3], type="l", main="Resources", xlab="Time", ylab="Concentration") matplot(r[,1], r[,4:12], type="l", main=paste(c("Spp 1","tolerance =",S1)), xlab="Time", ylab="Density") matplot(r[,1], r[,13:21], type="l", main=paste(c("Spp 2","tolerance =",S2)), xlab="Time", ylab="Density") #calculate R2 produced by sp1 R1N1<-sum(r[9991,4:12]*X*r[9991,2]/((f1*X)+r[9991,2])) #calculate R2 produced by sp2 R1N2<- sum(r[9991,13:21]*X*r[9991,2]/((f2*X)+r[9991,2])) #calculate R2 used by sp1 R2N1<- sum(r[9991,4:12]*(1-X)*r[9991,3]/((f1*(1-X))+r[9991,3])) #calculate R2 used by sp2 R2N2<- sum(r[9991,13:21]*(1-X)*r[9991,3]/((f2*(1-X))+r[9991,3])) # add to dataframe results[i,23]<-R1N1 results[i,24]<-R1N2 results[i,25]<-R2N1 results[i,26]<-R2N2 dev.off() print(i) } df<-data.frame(results[1:9,]) write.table(df, file='c:/folder1/results1.txt',row.names=FALSE) ################################# ## with allelopathic chemicals ## ################################# library(deSolve) bact.dynamics<-function(t,y,params) { with(as.list(params), { ## initial resource concentrations R1<-y[1] R2<-y[2] ## vector of density of each phenotype of the two species N1<-y[(length(R1)+length(R2)+1):(length(R1)+length(R2)+length(N1))] N2<-y[(length(R1)+length(R2)+length(N1)+1):(length(R1)+length(R2)+length(N1)+length(N2))] ## initial concentration of allelopathic chemicals A1<-y[(length(R1)+length(R2)+length(N1)+length(N2)+1)] A2<-y[(length(R1)+length(R2)+length(N1)+length(N2)+2)] #### if N falls below 0.01, call it 0 tmp<-which((N1)<0.01,arr.ind=T) tmp tmp2<-which((N1)>=0.01,arr.ind=T) tmp2 N1b<-array(NA, dim = length(N1)) N1b[tmp]<-(0) N1b[tmp2]<-(N1[tmp2]) N1<-N1b tmp<-which((N2)<0.01,arr.ind=T) tmp2<-which((N2)>=0.01,arr.ind=T) N2b<-array(NA, dim = length(N2)) N2b[tmp]<-(0) N2b[tmp2]<-(N2[tmp2]) N2<-N2b #### change in resource 1 is the sum of resources in through influx and removed by efflux minus those used by sp 1 and sp 2 #### vmax is a funtion of the trait X which has distribution traded-off with f (affinity for resource) and g (nATP produced in metabolism) from Gudelj et al 2008 #### sp1 and sp2 have different shaped tradeoff functions which are defined later dR1.dt<- dil.rate*(R1.food-R1) - sum((1/growth.yield.1)*N1*X*R1/((f1*X)+R1)) - sum((1/growth.yield.1)*N2*X*R1/((f2*X)+R1)) #### R2 is potentially produced and used by both species #### R2 is not brought in through influx but is lost in efflux #### growth on R2 traded-off against growth on R1 dR2.dt<- dil.rate*(-R2)+ sum((1/growth.yield.2)*factor.R2*N1*X*R1/((f1*X)+R1)) + sum((1/growth.yield.2)*factor.R2*N2*X*R1/((f2*X)+R1)) - sum((1/growth.yield.1)*N2*(1-X)*R2/((f2*(1-X))+R2)) - sum((1/growth.yield.1)*N1*(1-X)*R2/((f1*(1-X))+R2)) #### the growth of each spp will be limited by resource abundance so use resources/number of phenotypes in the later equations R1x<-c(R1/length(N1)) R2x<-c(R2/length(N2)) #### A1 is the amount of allelopathic chemical produced by species 1 # p is the probability of making A1 # this depends on the abundance of species 2 compared to the available resources P1<-sum(N2)/(sum(N1)+1) # if abundance of N2 is greater than total resources p=1. p cannot be greater than 1. p1<-ifelse(P1>1,1,P1) # c is probability of species 2 meeting an A1. This depends on how many A1s there are compared to individuals of species 2. c1<-A1*sum(N2)*eff #### change in A1 is the sum of A1 produced by sp1 minus those used against sp2 minus those removed by efflux dA1.dt<-eff*p1*sum(N1)-c1*A1*sum(N2)-A1*dil.rate # A1 cannot fall below 0 dA1.dt<- ifelse ((A1+dA1.dt)<0,(-A1),dA1.dt) #### A2 dynamics same as A1 P2<-sum(N1)/(sum(N2)+1) p2<-ifelse(P2>1,1,P2) c2<-A2*sum(N1)*eff dA2.dt<-eff*p2*sum(N2)-c2*A2*sum(N1)-A2*dil.rate # A2 cannot fall below 0 dA2.dt<- ifelse ((A2+dA2.dt)<0,(-A2),dA2.dt) #### The abundance of each phenotype is found by calculating growth on R1 plus that on R2 minus mortality caused by efflux, allelopatic chemicals and stress #### plus the number of individuals that are added to that phenotype by mutation of others (see Gudelj et al 2008) #### The susceptibility of a sp to allelopathic chemicals produced by others is traded off against the tolerance of that sp to stress # d is price to make an A dN1.dt<-array(NA, dim = length(N1)) dN2.dt<-array(NA, dim = length(N1)) # evolution # growth of N1 using R1 # growth of N1 using R2 # sp1 removed by efflux and stress (S1 is tolerance to stress) # price of producing A1s # mortality caused by A2 (S1 is the susceptibility of sp1 to A2) dN1a.dt<- mut.rate*(N1[2]-N1[1])+ c.conv*N1[1]*g1[1]*X[1]*(X[1]*R1x/((f1[1]*X[1])+R1x)) + c.conv*N1[1]*g1[1]*(1-X[1])*((1-X[1])*R2x/((f1[1]*(1-X[1]))+R2x)) - dil.rate*N1[1]/S1 - (eff*d*p1*N1[1])- (S1*c2*(A2)*N1[1]) dN1.dt[1]<-dN1a.dt for ( i in (2:(length(N1)-1))){ dN1i.dt<- mut.rate*(0.5*N1[i-1]+0.5*N1[i+1]-N1[i])+c.conv*N1[i]*g1[i]*X[i]*(X[i]*R1x/((f1[i]*X[i])+R1x)) + c.conv*N1[i]*g1[i]*(1-X[i])*((1-X[i])*R2x/((f1[i]*(1-X[i]))+R2x) ) -dil.rate*N1[i]/S1-(eff*d*p1*N1[i])-(S1*c2*(A2)*N1[1]) dN1.dt[i]<-dN1i.dt } dN1n.dt<- mut.rate*(N1[length(N1)-1]-N1[length(N1)]) +c.conv*N1[length(N1)]*g1[length(N1)]*X[length(N1)]*(X[length(N1)]*R1x/((f1[length(N1)]*X[length(N1)])+R1x)) + c.conv*N1[length(N1)]*g1[length(N1)]*(1-X[length(N1)])*((1-X[length(N1)])*R2x/((f1[length(N1)]*(1-X[length(N1)]))+R2x)) - dil.rate*N1[length(N1)]/S1 -(eff*d*p1*N1[length(N1)])-(S1*c2*(A2)*N1[length(N1)]) dN1.dt[length(N1)]<-dN1n.dt #### N2 dN2a.dt<- mut.rate*(N2[2]-N2[1])+ c.conv*N2[1]*g2[1]*(1-X[1])*((1-X[1])*R2x/((f2[1]*(1-X[1]))+R2x)) +c.conv*N2[1]*g2[1]*X[1]*(X[1]*R1x/((f2[1]*X[1])+R1x)) - dil.rate*N2[1]/S2-(eff*d*p2*N2[1])-(S2*c1*(A1)*N2[1]) dN2.dt[1]<-dN2a.dt for ( i in (2:(length(N2)-1))){ dN2i.dt<- mut.rate*(0.5*N2[i-1]+0.5*N2[i+1]-N2[i])+c.conv*N2[i]*g2[i]*(1-X[i])*((1-X[i])*R2x/((f2[i]*(1-X[i]))+R2x)) +c.conv*N2[i]*g2[i]*X[i]*(X[i]*R1x/((f2[i]*X[i])+R1x)) -dil.rate*N2[i]/S2 - (eff*d*p2*N2[i])-(S2*c1*(A1)*N2[i]) dN2.dt[i]<-dN2i.dt } dN2n.dt<- mut.rate*(N2[length(N2)-1]-N2[length(N2)]) +c.conv*N2[length(N2)]*g2[length(N2)]*(1-X[length(N2)])*((1-X[length(N2)])*R2x/((f2[length(N2)]*(1-X[length(N2)]))+R2x)) +c.conv*N2[length(N2)]*g2[length(N2)]*X[length(N2)]*(X[length(N2)]*R1x/((f2[length(N2)]*X[length(N2)])+R1x)) - dil.rate*N2[length(N2)]/S2-(eff*d*p2*N2[length(N2)])-(S2*c1*(A1)*N2[length(N2)]) dN2.dt[length(N2)]<-dN2n.dt #### when R falls below 0.01, call it 0 tmp<-which((R1)<0.01,arr.ind=T) tmp2<-which((R1)>=0.01,arr.ind=T) R1b<-array(NA, dim = length(R1)) R1b[tmp]<-(0) R1b[tmp2]<-(R1[tmp2]) R1<-R1b tmp<-which((R2)<0.01,arr.ind=T) tmp2<-which((R2)>=0.01,arr.ind=T) R2b<-array(NA, dim = length(R2)) R2b[tmp]<-(0) R2b[tmp2]<-(R2[tmp2]) R2<-R2b #### don't allow A to become negative tmp<-which((A1)<0.01,arr.ind=T) tmp2<-which((A1)>=0.01,arr.ind=T) A1b<-array(NA, dim = length(A1)) A1b[tmp]<-(0) A1b[tmp2]<-(A1[tmp2]) A1<-A1b tmp<-which((A2)<0.01,arr.ind=T) tmp2<-which((A2)>=0.01,arr.ind=T) A2b<-array(NA, dim = length(A2)) A2b[tmp]<-(0) A2b[tmp2]<-(A2[tmp2]) A2<-A2b return (list(c(dR1.dt, dR2.dt,dN1.dt,dN2.dt,dA1.dt,dA2.dt), NULL)) }) } #### specify parameters #dil.rate=0.05 ## rate of addition of resources and removal of resources, individuals and chemicals from the chemostat #eff=0.01 ## number of As produced by each N #R1.food=1000 ## concentration of resources in the influx #c.conv=0.75 ## conversion constant #factor.R2=1 ## how many R2s are produced per R1 used #mut.rate=0.05 ## mutation rate #f1=f1 ## vector describing the affinity of each phenotype of sp 1 for the resource #f2=f2 ## vector describing the affinity of each phenotype of sp 2 for the resource #g1=g1 ## vector describing the number of ATPs produced by each phenotype of sp 1 when using resource #g2=g2 ## vector describing the number of ATPs produced by each phenotype of sp 2 when using resource #X=X ## trait distribution #growth.yield.1=1 ## number of Ns produced by using R1 #growth.yield.2=1 ## number of Ns produced by using R2 #S1=S1 ## tolerance of N1 to stress (defined later). The inverse of S1 is the susceptibility of N1 to allelopathic chemical A2. #S2=S2 ## each phenotype has a given value of X which affects their vmax on R1 (vmax on R2 is a function of 1-X) ## species evolve by change in abundance of phenotypes X<-c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9) f1<-rep(0.1,9) f2<-rep(0.1,9) ## monod shaped tradeoff curves for g g1<-(X/(0.25+X)) g2<-1-(X/(0.75+X)) ## loop from here # data frame for results results<-matrix(NA,nrow=9,ncol=28) for (i in (1:9)){ # stress input values S1<-(i*0.1) S2<-0.9 params <- c(dil.rate=0.05,eff=0.01 ,R1.food=1000, c.conv=0.75, factor.R2=1, mut.rate=0.05,f1=f1,f2=f2,g1=g1,g2=g2,X=X,growth.yield.1=1,growth.yield.2=1,S1=S1,S2=S2,d=1) ## make initial populations of N1 and N2 v<-rnorm(9,2.5,sd=1) v<-ifelse(v<0,0,v) V<-sort(v) N1<-c(V[1],V[3],V[5],V[7],V[9],V[8],V[6],V[4],V[2]) v<-rnorm(9,2.5,sd=1) v<-ifelse(v<0,0,v) V<-sort(v) N2<-c(V[1],V[3],V[5],V[7],V[9],V[8],V[6],V[4],V[2]) # Initial values of resource, bacteria and allelopathic chemicals y <- c(r1=500,r2=0,n1=N1,n2=N2,A1=0,A2=0) # iterations times <- seq(1,1000,0.1) # Solving r <- lsoda(y, times, bact.dynamics, params,atol=0.0001,rtol=0.0001) # put results from final iteration into dataframe. Order: i*0.1, iteration number, abundance of R1, abundance of R2, abundance of sp1 when X=0.1, when X=0.2, when X=0.3 etc, abundance of sp2 when X=0.1, X=0.2 etc , abundance of A1, abundance of A2 results[i,1:24]<-c(i*0.1,as.vector(r[length(r[,1]),])) namefile<-c('c:/folder1/file1.pdf','c:/folder1/file2.pdf','c:/folder1/file3.pdf','c:/folder1/file4.pdf','c:/folder1/file5.pdf','c:/folder1/file6.pdf','c:/folder1/file7.pdf','c:/folder1/file8.pdf','c:/folder1/file9.pdf') pdf(file=namefile[i], height =6, width = 6) matplot(r[,1], r[,2:3], type="l", main="Resources", xlab="Time", ylab="Concentration") matplot(r[,1], r[,4:12], type="l", main=paste(c("Spp 1","tolerance =",S1)), xlab="Time", ylab="Density") matplot(r[,1], r[,13:21], type="l", main=paste(c("Spp 2","tolerance =",S2)), xlab="Time", ylab="Density") #calculate R2 produced by sp1 R1N1<-sum(r[9991,4:12]*X*r[9991,2]/((f1*X)+r[9991,2])) #calculate R2 produced by sp2 R1N2<- sum(r[9991,13:21]*X*r[9991,2]/((f2*X)+r[9991,2])) #calculate R2 used by sp1 R2N1<- sum(r[9991,4:12]*(1-X)*r[9991,3]/((f1*(1-X))+r[9991,3])) #calculate R2 used by sp2 R2N2<- sum(r[9991,13:21]*(1-X)*r[9991,3]/((f2*(1-X))+r[9991,3]) ) results[i,25]<-R1N1 results[i,26]<-R1N2 results[i,27]<-R2N1 results[i,28]<-R2N2 dev.off() print(i) } df<-data.frame(results[1:9,]) write.table(df, file='c:/folder1/results1.txt',row.names=FALSE)