Figure 10.14:

Tracking 10 particles undergoing different growth-rate dispersion mechanisms; intrinsic growth-rate dispersion (left), and growth-dependent dispersion (right).

Code for Figure 10.14

Text of the GNU GPL.

main.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
% comparing intrinsic growth rate dispersion to growth dependent dispersion
% jbr 10/5/04

nhalf = 5;
npart = 2*nhalf;

k = 10;

rand('seed',0)
randn('seed',0)

% particles have a distribution of growth rates
g0 = 1;
std = 0.1;
size0 = 10;
ntime = 15;
time = (0:ntime-1)';
x(:,1) = size0*ones(npart,1);
g = std*randn(npart,1) + g0;

for i = 1:ntime-1
    x(:,i+1) = x(:,i) + g;
end

% particles have same growth rate; choose which one grows at random
% particle increases size by increment delta

delta = 0.1;
nsim = 1500;
ipart = fix(rand(nsim,1)*npart)+1;
xran = size0*ones(npart,1);
simtim = (0:nsim-1)';
tau(1) = 0;

for i = 1:nsim-1
    tau(i + 1) = tau(i) - log(rand)/(k*npart);
    part = ipart(i);
    xran(:,i+1) = xran(:,i);
    xran(part,i+1) = xran(part,i) + delta;
end

%sort the particles by size; split into two groups, large and small,
%and grow them again

[size, label] = sort(x(:,ntime));

y1(:,1) = x(label(1:nhalf), ntime);
g1 = g(label(1:nhalf));
y2(:,1) = x(label(nhalf+1:npart), ntime);
g2 = g(label(nhalf+1:npart));

for i = 1:ntime-1
    y1(:,i+1) = y1(:,i) + g1;
    y2(:,i+1) = y2(:,i) + g2;
end

[size, label] = sort(xran(:,nsim));
y1ran(:,1) = xran(label(1:nhalf),nsim);
y2ran(:,1) = xran(label(nhalf+1:npart),nsim);
simtim = (0:(nsim-1)/2)';
tau1(1) = 0;
tau2(1) = 0;

ipart1 = fix(rand(nsim,1)*nhalf)+1;
ipart2 = fix(rand(nsim,1)*nhalf)+1;

for i = 1:(nsim - 1)/2
    tau1(i+1) = tau1(i) - log(rand)/(k*npart)*2;
    tau2(i+1) = tau2(i) - log(rand)/(k*npart)*2;
    part1 = ipart1(i);
    part2 = ipart2(i);
    y1ran(:,i+1) = y1ran(:,i);
    y2ran(:,i+1) = y2ran(:,i);
    y1ran(part1,i+1) = y1ran(part1,i) + delta;
    y2ran(part2,i+1) = y2ran(part2,i) + delta;
end

table1 = [time,x'];
table2 = [tau',xran'];

save dispmodel.dat table1 table2;

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
    subplot(1,2,1)
    plot (table1(:,1), table1(:,2:10));
    % TITLE dispmodel1

    subplot(1,2,2)
    plot (table2(:,1), table2(:,2:10));
    % TITLE dispmodel2
end % PLOTTING