1

Exercises Integrative Tumor Cell Biology

Analysis of tumor-growth models

B.W. Kooi, Dept. Theoretical Biology, FALW, Vrije Universiteit, Amsterdam

We will use program’s written in Matlab. Browse to Bbin Course Documents

and download the file ictbdynam.zip. In this file the scripts are organized intothree subdirectories. Download the fileictbdynam.zip in a folder in your home_dir.

Start MatlabwithStart/Programs/cursus/Matlab7.

In number of windows are opened. The Command Window is the main window in which you communicate with MATLAB

For each exercise, first change the current directory in the Current Directory field in the desktop toolbar. The files in the directory will appear in the Current Directory browser. For instance for analyzing the logistic growth model with a constant predation pressure, start the scriptrunlvp.mby double clicking in the listing of the Current Directory browser. Otherwise give the command directly after the prompt in the Command Window

> runlvp

Then give commands in the Command Window by answering questions.

Exercise1.Analyze the logistic growth model with constant predation pressure.

The differential equation is:

The first session will be:

> runlvp

K, carrying capacity (default value 1.0) is: 1

r growth rate (default value 1.0) is: 1

p predation preasure (default value p=0.2 (0<=p<0.5) is: 0

0<R(0) initial value R (default value 0.4) is: 0.4

starting point equilibrium finder:

R=0.999932

Func-count x f(x) Procedure

1 0.999932 6.81533e-005 initial

2 0.971649 0.0275468 search

3 1.02821 -0.0290102 search

Looking for a zero in the interval [0.97165, 1.0282]

4 0.9992 0.000799355 interpolation

5 1 -1.28264e-006 interpolation

6 1 1.02692e-009 interpolation

7 1 1.33227e-015 interpolation

8 1 0 interpolation

Zero found in the interval: [0.97165, 1.0282].

found equilibrium values:

R=1.000000

eigenvalue:

-1

do you want another run? y or n :

Typey if you want to run a new analysis with other parameter values, typeotherwise nand you get the prompt again.

The program calculates the equilibrium. In this example the equilibrium is R=1.000000 (the carrying capacity). Furthermore, the eigenvalue is calculated which indicates whether the equilibrium is stable or unstable (negative then stable, positive then unstable and if equal then zero labile).

In the figure below, the solution R(t) is plotted in the left panel. At equilibrium we have necessarily:

Therefore we plot in the right panel of the figure both the left- and the right-hand side of this equation as function of R. Their intersection point is an equilibrium or steady state.The vertical blue line gives the chosen initial value R(0).

The left figure gives the population size (number of tumor cells) R(t) as a function of the time t.

Repeat this session with various initial conditions or other parameter valuesas indicated below and answer the questions.

Do the same run but now with R(0)=0.2. Does the solution converge ultimately to the same final values as in the previous calculation?

The same now with p=0.2 andR(0)=0.2. To which value does the solution converge for large time? An equilibrium has been calculated. Which equilibrium is calculated and is it stable?

The same, now withR(0)=0.4. Is there convergence to the same value for R as in the pervious run? Is this equilibrium stable?

Is it possible to predict an unstable equilibrium by simulation?Does this hold also for other parameter values? For instance when two parameters are negative such that time changes from zero to minus infinity?

Now analyze the dynamical behavior with parameter valuep=0.25. Take again K=1, r=1 en p=0.25 and the initial conditionR(0)=0.6. To which point does the solution converge? Is there an equilibrium pointfound? Is this point stable?

The same, but now withp=0.25 andR(0)=0.4. What is the calculated equilibrium value? Is this equilibrium stable?

Exercise2.AnalyzetheGompertz model with constant predation pressure. The differential equation is:

This is the Gompertz equation whenp=0.

Change the current directory in the Current Directory field in the desktop toolbar to gompertz. The files will appear in the Current Directory browser. Start the scriptrungomp.m by double clicking in the listing of the Current Directory browser.

Remark: you always have to take K=1 and r0=1.

K, carrying capacity (default value 1.0) is: 1

r growth rate (default value 1.0) is: 1

p redation preasure (default value p=0.0 is: 0

0<R(0) initial value R (default value 0.1) is: 0.1

alpha =

2.3026

starting point equilibrium finder:

R=1.000000

R=0.100000

Func-count x f(x) Procedure

1 1 3.69431e-009 initial

2 0.971716 0.0641971 search

3 1.02828 -0.0660394 search

Looking for a zero in the interval [0.97172, 1.0283]

4 0.9996 0.000921095 interpolation

5 1 -4.91949e-007 interpolation

6 1 9.84489e-011 interpolation

7 1 0 interpolation

Zero found in the interval: [0.97172, 1.0283].

found equilibrium values:

R=1.000000

eigenvalue:

0

-2.3026

do you want another run? y or n :

The left panel gives the population size (number of tumor cells)R(t) as a function of the time t. At equilibrium we have necessarily:

Thereforewe plot in the right penal both the left- and the right-hand side of this equation as function ofR. Their intersection point is an equilibrium or steady state.The vertical line gives the chosen initial valueR(0).

Repeat the session with initialvalueR(0)=1.1. Do you agree with the result?

The same, but now with p=0.2 en R(0)=0.2. The horizontal line at the height of 0.2 indicates that there are two solutions. Obviously the pointR=0.86638 is stable. Do you think that the other equilibrium point is also stable?

The resulting figure is given here for p=0.2 en R(0)=0.2. Observe that the shape of the function

is different from the shape in the previous run with p=0en R(0)=0.1. Why?

Repeat de session with other initial values for example R(0)=0.001 and draw conclusions with respect to the stability of the other equilibrium point.

Exercise3. Analysis ofthe Goldbetter model of the cellcycle.

Minimum cascade model for a cell cycle oscillator reads:

where

Here C denotes the cyclin concentration, while M and X represent thefraction active cdc2 kinase and the fraction inactive cdc2 kinase.

Change the current directory in the Current Directory field in the desktop toolbar to goldcell.

> rungold

K_1, K_2 (default value 0.005), is:

K_1, K_2=0.005000

K_3, K_4 (default value 0.005), is:

K_3, K_4=0.005000

0<C(0) (default value 0.01) is:

C(0)=0.010000

0<M(0)<1 (default value 0.01) is:

M(0)=0.010000

0<X(0)<1 (default value 0.01) is:

X(0)=0.010000

starting point equilibrium finder:

C=0.507973, M=0.158681, X=0.002286

k= 2, C=1.026302, M=2.945228, X=0.061180

k= 3, C=0.199583, M=0.484509, X=0.094752

k= 4, C=0.409295, M=0.474542, X=0.082940

k= 5, C=0.490931, M=0.474332, X=0.083519

k= 6, C=0.498914, M=0.474245, X=0.083251

k= 7, C=0.498979, M=0.474244, X=0.083249

k= 8, C=0.498979, M=0.474244, X=0.083249

found equilibrium values:

C=0.498979, M=0.474244, X=0.083249

eigenvectors:

eigenvalues:

0.2300 + 0.6032i

0.2300 - 0.6032i

-0.8546

do you want another run? y or n :

You see four figures. The upper two panelsare phase-plain plots in whichM andX versusC are plotted. The figure panelleft-bottom gives the three state variablesC,M,X as function of timet. The figure panel right-bottom is a two-parameter diagram. The vertical and the horizontal lines indicate the values of the parameters K1 =K2 and the parameters K3 = K4.In the region right-above the solid curve the system possesses a stable equilibrium. Their intersection is the chosen point in this diagram. In our case the point is left-below the solid curve. In this region the equilibria are unstable and the solution oscillates, see also left-bottompanel.

Twoof the calculatedeigenvalues are complex, 0.23001 + 0.60320i, whereby thefirst part is the real part and the second part is the imaginary part where i=√-1.These two eigenvalues possess positive real part, and therefore the equilibrium is unstable.

Repeat the calculation, now withK1=K2 =1 en K3 = K4 =1. Is the equilibrium stable? Are the real parts of all three eigenvalues negative? Is in the figure right-bottom the intersection of the horizontal and thevertical line in the region where the equilibrium is stable?