References:

1. http://www.math.uic.edu/~tier/Math419/Tumor/Notes/gompertz.html

Abstract. We provide a numerical FORTRAN code to solve the system of DE known as Gompertz equations.

The Gompertz equation can be used to describe the growth  of numerous tumors.

The volume of a tumor may grow at a rate r ~ 1/time

  d V(t) /dt = r(t) V(t)    , with specified initial volume    V(0)= V₀ ~ length³                       .     ( 1)

However the rate r may also decay in time like

 d r(t) /dt = - α r(t )       α ~ 1/time  with      IC   r(0) = β ~ 1/time       .                          (2)

The analytical solution is , where α is assumed to be constant ,

                                    V(t) = V₀ exp {(β /α) ( 1- exp(-αt) ) }             .  ( 3)

There two time scales in the problem

 t_scale1 = 1/α    , t_scale2 = 1/β .  Equations (1) and (2) have solutions

V(t_n) = V(t_n-1 ) +  (∆t) r (t_n-1 ) V(t_n-1 )                                          ,      ( 4)

r(t_n) = r(t_n-1 )  - α  (∆t)  r(t_n-1 )                                               .            (5)

The choice of ∆t must be smaller than the smallest of   t_scale1  and  t_scale2 .

Example :

Let  α  =1/.01 days^-1 , β = .04 days^-1  and  V₀ = 5.0 cm³ .Find V(t) , solving  the set of eqs (4) -(5) .

The time scales are  t_scale1 = 100 days   ,  t_scale2 = 25 days .

The results show a  growth from V₀ = 5 cm³ (radius ≈ 1.1 cm )   to  V = 253 cm³ ( radius ≈ 3.9 cm ) in 399 days.

The numerical solution matches very well the analytical calculation.

Fig 1. Tumor growth according to Gompertz equation.

FORTRAN code

c The Gomperz equation
c Ref http://www.math.uic.edu/~tier/Math419/Tumor/Notes/gompertz.html
c dV/dt = r(t) * V(t) ; dr/dt = alpha*r(t) ;IC V(0)-V0 ; r(0)= beta
c alpha ~ 1/time , beta ~ 1/time
data alpha , beta , vzero/.01 , .04 , 5. /
vgom(t)= vzero*exp( (beta/alpha)*(1.-exp(-alpha*t) ) )
tscale1=amin1(1./alpha,1./beta)
tscale2=amax1(1./alpha,1./beta)
dt=tscale1/1000.
tfinal=4.*tscale2
nstep=int(tfinal/dt)
kp=int(float(nstep)/50.)
kount=kp
r0=beta
v0=vzero
print 100 ,0., v0 ,vgom(0.)
do 10 i=1,nstep
t=dt*float(i)
v1=v0+dt*(r0*v0)
r1=r0-dt*alpha*r0
if(i.eq.kount)then
print 100 ,t, v1 ,vgom(t)
kount=kount+kp
endif
v0=v1
r0=r1
10 continue
100 format('t,v,Vgompertz=',3(3x,e10.3))
stop
end