The following code constituted the fortran program, comp10.
c This simulates the initial infection
c by HBV, spread to the entire liver and growth of copy number to
c assigned limits. During viral growth some cells may become
c integrated. The program continues immune cell killing and so
c abets cell division and subsequent viral growth. This growth
c increases the fraction of the liver with viral integrants. Only
c one integrant per cell will be considered. March 25, 2008
c July 17, 2008. Added pzero1,pzero2, early and late chances
c of losing ccc DNA during mitosis, (i.e. before/after nday1)
implicit none
logical prt,growing
character ans,ans2,ans3,ans2b,fracfix
real repfrac,fracgrow,fixgrow
real reprate,pdinf0,pdinf1,pzero1,pzero2,eps
real pzero,comp,livereq,avcopies,efficiency,tcomp,pint
real finf,fboth,fintninf,finfnint,fnone,lambda,pdinf,pdu
integer ntarget,nncells,maxncells
integer ncells,ndays,nd,lower,upper,i
integer c(999999),copydist(0:60),nday1
integer n1,nkillcells,ijk,nn1,nn2,nn3,nn4,sumt
common ijk,ans3
eps=1.e-6
type*,'Want a new set of random numbers? (y/n)'
accept 61,ans3
61 format(a1)
call getparams(lower,upper,lambda,ans,pdinf1,pdu,ndays,
* finf,fboth,fintninf,ncells,ans2,pzero1,efficiency,
* pint,nday1,reprate,pdinf0,pzero2,repfrac,fracfix,
* fracgrow,fixgrow)
maxncells=ncells
finfnint=finf-fboth
fnone=1.-fboth-fintninf-finfnint
growing=.false.
pdinf=pdinf0
call initialize(ncells,c,ans,lower,upper,lambda,
* fboth,fintninf,finfnint,fnone)
nkillcells=pdinf1*finf*maxncells ! for use with fixed kill number
if(ans2.eq.'y')then
type*,'number of cells killed per day',nkillcells
c pause
endif
call account(c,ncells,nn1,nn2,nn3,nn4,sumt)
c open(unit=1,file='comp10.txt',status='unknown')
open(unit=1,file='comp10.csv',status='unknown')
nd=0
livereq=0.
prt=.true.
call clones(c,ncells,efficiency,comp,prt,nd,tcomp)
call copies(c,ncells,avcopies,copydist)
write(1,59) nd,comp,
* avcopies,livereq/ncells,nn1,nn2,nn3,nn4,sumt,tcomp
c 59 format(x,i7,3x,3(f15.10,3x),4(i7,x),i7,x,f10.6)
59 format(x,i7,',',3(f15.10,','),5(i7,','),f10.6)
c Bills version is comma delimited
ans2b='n'
do nd=1,ndays
pzero=pzero1
if(nd.gt.nday1)then
pdinf=pdinf1
pzero=pzero2
if(ans2.eq.'y')ans2b='y'
endif
type*,'day: ',nd
if(ans2b.ne.'y')then
call killcells(ncells,n1,c,pdinf,pdu) ! n1 cells left
type*,'nd,ncells,n1,ncells-n1',nd,ncells,n1,ncells-n1
else
type*,'nd,ncells,n1,nkillcells',nd,ncells,n1,nkillcells
call fixedkill(ncells,n1,nkillcells,c,pdu)
endif
livereq=livereq+ncells-n1
ncells=n1
nncells=repfrac*maxncells
if(n1.le.nncells.or.growing)then
growing=.true.
if(fracfix.eq.'f')then
ntarget=n1*(1.+fracgrow) ! grows a fraction of its
else ! current size
ntarget=n1+fixgrow*maxncells ! grows a fraction of
endif ! the full liver
ntarget=min(1+ntarget,maxncells)
c if(ntarget.eq.maxncells)growing=.false. !only one growth cycle
call replicate(c,ntarget,n1,pzero) !above line is commented out
ncells=ntarget
c call replicate(c,ncells,n1,pzero)
endif
if(nd.le.nday1)then
call upregulate(c,ncells,pint,reprate)
if(reprate.gt.eps)then
call infect(c,ncells) ! puts one virus in any uninfected cell
endif
endif
call account(c,ncells,nn1,nn2,nn3,nn4,sumt)
prt=.false.
if(nd.eq.ndays.or.nd.eq.nday1.or.nd.eq.1)prt=.true.
call clone(sc,ncells,efficiency,comp,prt,nd,tcomp)
call copies(c,ncells,avcopies,copydist)
write(1,59) nd,comp,
* avcopies,livereq/maxncells,nn1,nn2,nn3,nn4,sumt,tcomp
c open(unit=33,file='cpydist10.txt',status='unknown')
open(unit=33,file='cpydist10.csv',status='unknown')
if(nd.eq.1.or.mod(nd,10).eq.0.or.nd.eq.nday1.or.
* nd.eq.ndays)then
type*,' '
if(nd.eq.1)write(33,*)'day',',','copies',',','cells'
type*,'distribution (no zeros) of copy number on day:',nd
do i=0,60
if(copydist(i).gt.0)then
type*,nd,i,copydist(i)
write(33,50)nd,i,copydist(i)
50 format(x,i7,',',i7,',',i7)
c 50 format(x,i7,'',i7,'',i7)
endif
enddo
type*,' '
endif
enddo
close(unit=1)
type*,' '
type*,'final distribution of copy number'
type*,' day copies cells'
do i=0,60
if(copydist(i).gt.0)then
type*,ndays,i,copydist(i)
endif
60 format(x,i4,5x,i6,5x,i10)
enddo
type*,' '
open(unit=13,file='seed.dat',status='unknown')
write(13,*)ijk
close(unit=13)
type*,' '
type*,'output file: comp10.csv, contains the following:'
type*,'day, complexity, average cccDNA copies, liver
* equivalent cell deaths, int&inf,int¬-inf,not-int&inf,
* not-int¬-inf, total cells, true complexity'
c n1: integrated and infected
c n2: integrated and not-infected
c n3: not integrated and infected
c n4: not integrated and not-infected
type*,' '
type*,'output file: cpydist10.csv contains day,copy number,
* number of cells with this copy number. (every 10 days)
* days are in serial order, ie. the whole file is two columns.'
type*,' '
type*,'output file: clonesize10.csv, contains size and number of
* clones of this size'
type*,' '
type*,'input and output file: seed.dat -> random number seed'
type*,' '
close(unit=33)
pause 'hit enter to terminate program'
end
subroutine copies(c,ncells,avcopies,copydist)
implicit none
real avcopies
integer copydist(0:60),i,target
integer c(1),ncells,n,id,cnum,copynum
do i=0,60
copydist(i)=0
enddo
avcopies=0.
do n=1,ncells
cnum=copynum(c,n,id,target)
avcopies=avcopies+cnum
copydist(cnum)=copydist(cnum)+1
enddo
avcopies=avcopies/ncells
return
end
function b(cnum) ! random division of ccc DNA copies
implicit none ! to one daughter cell
integer b,cnum,i,idum
real rn
b=0
do i=1,cnum
if(rn(idum).le.0.5)b=b+1
enddo
return
end
subroutine replicate(c,ncells,n1,pzero)
implicit none
integer cnum,copynum,id,cn1,cn2,b
integer idum,need,i,j,temp,target
integer c(1),n1,k1,ncells
real pzero,rn
need=ncells-n1
k1=n1
c type*,'k1 at start of replicate',k1
do i=1,need
j=1+k1*rn(idum)
cnum=copynum(c,j,id,target)
cnum=nint(cnum*(1.-pzero)) ! loss of copies during mitosis
cn1=b(cnum)
cn2=cnum-cn1 ! cn1 and cn2 will both be zero if c was not
temp=c(k1) ! an infected cell
c c(k1)=id+1000000*cn1
c type*,'*before newcopynum: k1,id,cn1,target', k1,id,cn1,target
call newcopynum(c,k1,id,cn1,target)
c type*,'* after newcopynum: k1,id,cn1,target', k1,id,cn1,target
c subroutine newcopynum(c,n,id,cnum,target)
c(j)=temp ! allows only one replication per cell per day
c c(n1+i)=id+1000000*cn2
call newcopynum(c,n1+i,id,cn2,target)
k1=k1-1
enddo
return
end
subroutine fixedkill(ncells,n1,killcells,c,pdu) ! n1 cells left
implicit none
integer j,nkill,nkinf,nku,nk,cnum,id,target,copynum
integer c(1),idum,n1,n,ncells,killcells,killset(999999)
real pdu,rn
n1=0
j=1
do n=1,ncells
cnum=copynum(c,n,id,target)
if(cnum.ge.1)then ! collect infected cells into killset
killset(j)=c(n)
j=j+1
endif
enddo
nkill=j-1
nk=nkill
if(killcells.le.nkill)then
do n=1,killcells
j=1+nk*rn(idum)
killset(j)=killset(nk)
killset(nk)=-1
nk=nk-1
enddo ! killset now contains all infected cells that are alive
nkinf=nkill-killcells ! number of live cells in killset
else
nkinf=0
endif
j=1
do n=1,ncells
cnum=copynum(c,n,id,target)
if(cnum.eq.0)then ! kill uninfected cells
if(rn(idum).gt.pdu)then
c(j)=c(n)
j=j+1
endif
endif
enddo
nku=j-1 ! number of live uninfected cells, now compressed 1 to nku
do n=nku+1,nku+nkinf
c(n)=killset(n-nku)
enddo
do n=nku+nkinf+1,ncells
c(n)=-1
enddo
n1=nkinf+nku
return
end
subroutine killcells(ncells,n1,c,pdinf,pdu)
implicit none
integer c(1),idum,n1,n,ncells,nlast,cnum,copynum,target,id
real pdinf,pdu,rn
n1=0
do n=1,ncells
cnum=copynum(c,n,id,target)
if(cnum.ge.1)then
if(rn(idum).le.pdinf)then
c(n)=-1
endif
else
if(rn(idum).le.pdu)then
c(n)=-1
endif
endif
enddo
nlast=ncells
n1=0
do n=ncells,1,-1
if(c(n).eq.-1)then
n1=n1+1
c(n)=c(nlast)
c(nlast)=-1
nlast=nlast-1
endif
enddo
n1=ncells-n1
return
end
subroutine account(c,ncells,n1,n2,n3,n4,sumt)
implicit none
integer sumt,n1,n2,n3,n4,cnum,id
integer c(1),i,ncells,copynum,target
c n1: integrated and infected
c n2: integrated and not-infected
c n3: not integrated and infected
c n4: not integrated and not-infected
n1=0
n2=0
n3=0
n4=0
do i=1,ncells
cnum=copynum(c,i,id,target)
if(cnum.gt.0.and.id.gt.0)n1=n1+1 ! infected, integrated
if(cnum.eq.0.and.id.gt.0)n2=n2+1 ! not infected, integrated
if(cnum.gt.0.and.id.eq.0)n3=n3+1 ! infected, not integrated
if(cnum.eq.0.and.id.eq.0)n4=n4+1 ! uninfected, integrated
enddo
c type*,'integrated and infected cells ',n1
c type*,'integrated not infected cells ',n2
c type*,'not integrated and infected cells ',n3
c type*,'not integrated not infected cells ',n4
sumt=n1+n2+n3+n4
c type*,'overall liver size',sumt
return
end
subroutine initialize(ncells,c,ans,lower,upper,lambda,
* fboth,fintninf,finfnint,fnone)
implicit none
integer n1,n2,n3,n4,idum
integer ncells,c(1),lower,upper,n,cnum,npois
real lambda,fboth,fintninf,finfnint,fnone,rn
character ans
n1=nint(fboth*ncells) ! both integrated and infected
do n=1,n1 ! both integrated and infected ! 0 (zero) is the first ID.
if(ans.eq.'u')then ! u for uniform distribution
cnum=lower+(1+upper-lower)*rn(idum)
else
cnum=npois(lambda,lower,upper) ! ans = 'p'
endif
call newcopynum(c,n,n,1,cnum)
c c(n)=n+cnum*1000000
enddo
c subroutine newcopynum(c,n,id,cnum,target)
n2=nint(fintninf*ncells)
do n=n1+1,n1+n2 ! top down
c c(n)=n ! integrated but not infected (zero ccc DNA copies)
if(ans.eq.'u')then ! u for uniform distribution
cnum=lower+(1+upper-lower)*rn(idum)
else
cnum=npois(lambda,lower,upper) ! ans = 'p'
endif
call newcopynum(c,n,n,0,cnum) ! has potential to be infected
enddo
n3=nint(finfnint*ncells) ! infected but not integrated
do n=n1+n2+1,n1+n2+n3
if(ans.eq.'u')then
cnum=lower+(upper-lower)*rn(idum)
else
cnum=npois(lambda,lower,upper) ! ans = 'p'
endif
call newcopynum(c,n,0,1,cnum)
c c(n)=cnum*1000000 ! cell id is zero
enddo
n4=ncells-n1-n2-n3
do n=n1+n2+n3+1,ncells
if(ans.eq.'u')then
cnum=lower+(upper-lower)*rn(idum)
else
cnum=npois(lambda,lower,upper) ! ans = 'p'
endif
call newcopynum(c,n,0,0,cnum) ! has potential to be infected
c c(n)=0
enddo
return
end
real function rn(idum)
implicit none
character ans3 ! read seed or start from scratch (y/n)
real rug(101),ran_,xxx
integer init,ijk,idum,i,j
common ijk,ans3
data init/0/
if(init.eq.0)then
init=1
if(ans3.eq.'y')then
open(unit=13,file='seed.dat',status='old')
read(13,*)ijk
close(unit=13)
type*,'read old seed and will continue from there'
endif
do i=1,3389
xxx=ran_(ijk)
enddo
do i=1,101
rug(i)=ran_(ijk)
enddo
endif
j=1+100.d0*ran_(ijk)
rn=rug(j)
rug(j)=ran_(ijk)
return
end
integer function npois(lambda,lower,upper)
implicit none
real lambda,a,p,rn
integer idum,init,lower,upper
data init/0/
if(init.eq.0)then
init=1
a=exp(-lambda)
endif
1 p=rn(idum)
npois=0
dowhile(p.gt.a)
p=p*rn(idum)
npois=npois+1
enddo
if(npois.lt.lower.or.npois.gt.upper)go to 1
return
end
subroutine getparams(lower,upper,lambda,ans,pdinf,pdu,ndays,
* finf,fboth,fintninf,ncells,ans2,pzero1,efficiency,
* pint,nday1,reprate,pdinf0,pzero2,repfrac,fracfix,
* fracgrow,fixgrow)
implicit none
integer upper,lower,ndays,ncells,k,nday1
real pint,reprate,pdinf0,repfrac,fracgrow,fixgrow
real lambda,pdinf,pdu,finf,fboth,fintninf,pzero1,pzero2,efficiency
character ans,ans2,fracfix
save
open(unit=2,file='comp10.in',status='old')
read(2,*)lower
read(2,*)upper
read(2,*)lambda
read(2,51)ans ! u/p
read(2,*)pdinf
read(2,*)pdu
read(2,*)ndays ! 7
read(2,*)finf
read(2,*)fboth ! 9
read(2,*)fintninf
read(2,*)ncells
read(2,51)ans2
read(2,*)pzero1,pzero2
read(2,*)efficiency
read(2,*)pint
read(2,*)nday1
read(2,*)reprate
read(2,*)pdinf0
read(2,*)repfrac
read(2,*)fracfix
read(2,*)fracgrow
read(2,*)fixgrow
ncells=min(ncells,999999)
close(unit=2)
51 format(a1)
53 format(' 4. Uniform or Poisson distribution (u/p): ',a1)
1 type*,' '
type*,' 1. lower limit of initial CCCdna copies per cell: ',lower
type*,' 2. upper limit of initial CCCdna copies per cell: ',upper
type*,' 3. mean of Poisson distribution of CCCdna copies: ',lambda
type 53,ans ! '4. Uniform or Poisson distribution (u/p): '
type*,' 5. Fraction of infected cells dying per day: ', pdinf
type*,' 6. Fraction of UN-infected cells dying per day: ', pdu
type*,' 7. Number of days (cycles) to be computed: ',ndays
type*,' 8. Fraction of cells initially infected: ',finf
type*,' 9. Initial fraction of infected cells that are integrated:
* ',fboth
type*,'10. Initial fraction of UNinfected that are integrated: ',
* fintninf
type*,'11. number of cells in simulated liver (up to a million)'
* ,ncells
type 54,ans2 ! Fixed number of cells killed per day? (y/n)
54 format(x,'12. Fixed number of cells killed per day? (y/n): ',a1)
type*,'13. Frctns DNA lost on rep (early/late): ',pzero1,pzero2
type*,'14. Efficiency of integrant detection: ',efficiency
type*,'15. Chance of integration per viral replication event: ',
* pint
type*,'16. Number of days of UNinhibited viral growth',nday1
type*,'17. Viral rep rate per copy per day: ',reprate
type*,'18. Initial fraction infected cells dying per day',pdinf0
type*,'19. Liver fraction when replication kicks in',repfrac
type*,'20. Fractional or fixed growth per day?(f/x) ',fracfix
type*,'21. Fractional growth per day (if f)',fracgrow
type*,'22. Fixed growth per day (if x)',fixgrow
type*,' '
type*,'enter 0 if done, otherwise choose a number in [1,...,22]'
accept*,k
select case(k)
case(0)
open(unit=2,file='comp10.in',status='unknown')
write(2,*)lower
write(2,*)upper
write(2,*)lambda
write(2,51)ans
write(2,*)pdinf
write(2,*)pdu
write(2,*)ndays
write(2,*)finf
write(2,*)fboth
write(2,*)fintninf
write(2,*)ncells
write(2,51)ans2
write(2,*)pzero1,pzero2
write(2,*)efficiency
write(2,*)pint
write(2,*)nday1
write(2,*)reprate
write(2,*)pdinf0
write(2,*)repfrac
write(2,*)fracfix
write(2,*)fracgrow
write(2,*)fixgrow
close(unit=2)
open(unit=2,file='legend.out',status='unknown')
write(2,*)'lower limit: ',lower
write(2,*)'upper limit: ',upper
write(2,*)'mean copies: ',lambda
write(2,553)ans ! '4. Uniform or Poisson distribution (u/p): '
553 format(x,'Unif or Pois: ',a1)
write(2,*)'Kd-infected: ', pdinf
write(2,*)'Kd-uninfected: ', pdu
write(2,*)'days: ',ndays
write(2,*)'fraction infctd:',finf
write(2,*)'infctd intgrtd: ',fboth
write(2,*)'UNnfctd intgrtd: ',fintninf
write(2,*)'liver size: ',ncells
write(2,554)ans2 ! Fixed number of cells killed per day? (y/n)
554 format(x,'Fixed cell kill? : ',a1)
write(2,*)'DNA lost in replctn: ',pzero1,pzero2
write(2,*)'Efficiency of integrant detection: ',efficiency
write(2,*)'Chance of integration per replication: ',pint
write(2,*)'days of uninhibited viral growth: ',nday1
write(2,*)'viral rep rate per copy per day: ',reprate
write(2,*)'initial infected cell kill rate ',pdinf0
write(2,*)'Liver fraction when replication kicks in',repfrac
write(2,555)fracfix
555 format(x,'Frction or fixed liver growth per day:',a1)
write(2,*)'Fractional growth per day (if f)',fracgrow
write(2,*)'Fixed growth per day (if x)',fixgrow
close(unit=2)
return
case(1)
type*,'1. Lower limit of copy number initial distribution:'
accept *,lower
case(2)
type*,'2. Upper limit of copy number initial distribution: '
accept*,upper
case(3)
type*,'3. enter average copy number'
accept*,lambda
case(4)
type*,'4. enter u or p for uniform or poisson distribution'
accept 51,ans
case(5)
type*,'5. kd for infected cells: (fraction killed per day)'
accept*,pdinf
case(6)
type*,'6. kd for uninfected cells: '
accept*,pdu
case(7)
type*,'7. enter number of days: '
accept*,ndays
case(8)
type*,'8. enter fraction of cells initially infected: '
accept*,finf
case(9)
type*,'9. enter fraction of initially infected cells that
* are also integrated: '
accept *,fboth
case(10)
type*,'10. enter fraction UNinfected that are integrated: '
accept*,fintninf
case(11)
type*,'11. enter number of cells in simulated
*liver (at most 999,999)'
accept*,ncells
if(ncells.gt.999999)ncells=999999
case(12)
type*,'12. fixed number of cells killed per day? (y/n) '
accept 51,ans2
case(13)
type*,'13. enter frctns ccc DNA lost on rep (early/late): '
accept*,pzero1,pzero2
case(14)
type*,'14. enter efficiency of integrant detection: '
accept*,efficiency
case(15)
type*,'15. Enter chance of integration per replication: '
accept*,pint
case(16)
type*,'16. Enter days of uninhibited viral growth: '
accept*,nday1
case(17)
type*,'17. Enter viral replication rate per day per copy: '
accept*,reprate
case(18)
type*,'18. Enter initial kill rate of infected cells: '
accept*,pdinf0
case(19)
type*,'19. Enter Liver fraction when replication kicks in:'
accept*,repfrac
case(20)
type*,'20. Fractional or fixed liver growth per day? (f/x)'
accept 556,fracfix
556 format(a1)
case(21)
type*,'21. Enter Fractional growth per day (if f)'
accept*,fracgrow
case(22)
type*,'22. Enter Fixed growth per day (if x)'
accept*,fixgrow
case default
type*,'illegal entry, try again'
end select
go to 1
end
subroutine clones(c,ncells,efficiency,comp,prt,nd,tcomp)
c distribution of clone sizes
implicit none
logical prt ! to write out copynum files
integer i,clonesize(0:100000),k,names(0:999999),idum,nd,target
integer c(1),ncells,n,id,cnum,copynum,numclones,init
integer tclonesize(0:100000),tcnum,tnumclones,tnames(0:999999),kt
real efficiency,rn,comp,tcomp
data init/0/
if(init.eq.0)then
c open(unit=17,file='clonesize10.txt',status='unknown')
open(unit=17,file='clonesize10.csv',status='unknown')
init=1
endif
do i=1,100000
clonesize(i)=0
tclonesize(i)=0
enddo
do i=1,999999
names(i)=0
tnames(i)=0
enddo
do n=1,ncells
cnum=copynum(c,n,id,target)
tnames(id)=tnames(id)+1
if(rn(idum).le.efficiency)names(id)=names(id)+1
enddo
do id=1,ncells
kt=min(100000,tnames(id))
k=names(id)
k=min(k,100000)
tclonesize(kt)=tclonesize(kt)+1
clonesize(k)=clonesize(k)+1 ! counts number of clones of size k
enddo
if(prt)then
type*,'size and number of clones of this size'
type*,'size observed number true number'
endif
cnum=0
tcnum=0
if(prt)then
c open(unit=17,file='clonesize10.txt',status='unknown')
open(unit=17,file='clonesize10.csv',status='unknown')
write(17,*)' '
write(17,*)'day:',nd
write(17,48)'size',',','observed number',',','true number'
endif
48 format(x,a4,a1,a15,a1,a11)
49 format(x,i7,',',i7,',',i7)
numclones=0
tnumclones=0
do i=1,100000
if(tclonesize(i).gt.0)then
numclones=numclones+clonesize(i)
tnumclones=tnumclones+tclonesize(i)
if(prt)then
type*,i,clonesize(i),tclonesize(i)
write(17,49)i,clonesize(i),tclonesize(i)
endif
cnum=cnum+i*clonesize(i)
tcnum=tcnum+i*tclonesize(i)
endif
enddo
c close(unit=17)
if(prt)then
type*,'(obs/true) distinct clones',numclones,tnumclones
type*,'(obs/true) integrated cells ',cnum,tcnum
endif
comp=numclones
if(cnum.eq.0)then
comp=0.
else
comp=comp/cnum
endif
tcomp=tnumclones
if(tcnum.eq.0)then
tcomp=0.
else
tcomp=tcomp/tcnum
endif
if(prt)then
type*,'day',nd,' observed complexity',comp,' true complexity',
* tcomp,' (zero if undefined)'
endif
return
end
subroutine newcopynum(c,n,id,cnum,target)
implicit none ! submit n,id,cnum and target
integer itarget,cnum,knum,j ! load cell c(n)
integer invtarget(0:60,0:60),init
integer c(1),n,id,target
data init/0/
if(init.eq.0)then
j=1
do itarget=0,60
do knum=0,itarget
invtarget(knum,itarget)=j
j=j+1
enddo
enddo
init=1
endif
c(n)=id+1000000*invtarget(cnum,target)
return
end
function copynum(c,n,id,target)
implicit none ! submit n, get id,target and copynumber
integer itarget,cnum,j,knum
integer ntarget(0:2140,2),init
integer c(1),copynum,n,id,k,target
data init/0/
if(init.eq.0)then
j=1
do itarget=0,60
do knum=0,itarget
ntarget(j,2)=itarget
ntarget(j,1)=knum
j=j+1
enddo
enddo
init=1
endif
k=c(n)
cnum=k/1000000
id=k-cnum*1000000
if(cnum.gt.2140.or.cnum.lt.0)type*,'cnum > 2140 or < 0',cnum
copynum=ntarget(cnum,1)
target=ntarget(cnum,2)
return
end
subroutine upregulate(c,ncells,pint,reprate)
implicit none
integer idum,maxid,id,newpois,newcopies,init
integer c(1),ncells,copynum,cnum,n,target
real pint,reprate,rn
data init/0/
if(init.eq.0)then
maxid=0
do n=1,ncells
cnum=copynum(c,n,id,target)
maxid=max(id,maxid)
enddo
type*,'maxid on first call',maxid
init=1
endif
do n=1,ncells
cnum=copynum(c,n,id,target)
if(cnum.lt.target.and.cnum.gt.0)then ! upregulate
c type*,'cnum,newcopies,target -',cnum,newcopies,target
newcopies=newpois(cnum,reprate,target-cnum)
c type*,'cnum,newcopies,target +',cnum,newcopies,target
if(id.eq.0)then ! cell not integrated
if(rn(idum).le.newcopies*pint)then
maxid=maxid+1
id=maxid
if(maxid.gt.999999)stop 'id numbering exceeded'
endif
endif
call newcopynum(c,n,id,cnum+newcopies,target)
endif
enddo
return
end
integer function newpois(cnum,reprate,upper)
implicit none ! cum(i,iupper,icopies)
c logical done
c real fac(0:60),sum(60,60),xxx
real lambda,a,p,rn,reprate
c real mean,cum(0:60,1:60,1:60)
integer idum,upper,cnum
c integer idum,upper,n,init,icopies,iupper,i,j,cnum,mupper
c data init/0/
lambda=cnum*reprate
a=exp(-lambda)
1 p=rn(idum)
newpois=0
dowhile(p.gt.a)
p=p*rn(idum)
newpois=newpois+1
enddo
if(newpois.gt.upper)newpois=upper
return
c The version below may be useful later.
c if(init.eq.0)then
c init=1
c fac(0)=0.
c do i=1,60
c fac(i)=fac(i-1)+log(real(i))
c enddo
c do icopies=1,60
c mean=log(reprate)+log(real(icopies))
c do iupper=1,60
c xxx=0.
c do i=0,60
c cum(i,iupper,icopies)=exp(-mean+j*mean-fac(j))
c xxx=xxx+cum(i,iupper,icopies)
c enddo
c sum(icopies,iupper)=xxx
c do i=0,10
c cum(i,iupper,icopies)=cum(i,iupper,icopies)/xxx
c enddo
c enddo
c enddo
c endif
c xxx=rn(idum)
c i=0
c done=.false.
c do while(.not.done.and.i.le.upper)
c if(xxx.le.cum(i,upper,cnum))then
c done=.true.
c newpois=i
c else
c i=i+1
c endif
c enddo
c return
end
subroutine infect(c,ncells)
implicit none
integer c(1),ncells,n,copynum,target,cnum,id
do n=1,ncells
cnum=copynum(c,n,id,target)
if(cnum.eq.0)call newcopynum(c,n,id,1,target)
enddo
return
end