TVE-F 17 013 juni
Examensarbete 15 hp Juni 2017
Mathematical modeling of interactions between colonic crypts
Nils Erlanson
Martin Edin
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress:
Box 536 751 21 Uppsala Telefon:
018 – 471 30 03 Telefax:
018 – 471 30 00 Hemsida:
http://www.teknat.uu.se/student
Abstract
Mathematical modeling of interactions between colonic crypts
Nils Erlanson, Martin Edin
Mathematical modeling of colonic crypts has been done for almost half a century, partly as corectal cancer is widely suspected to originate in the crypts. This paper will investigate how cells proliferate in colonic crypts, and the interaction with adjacent crypts by mathematical modeling. The model uses a pressure gradient to simulate the proliferation. The pressure gradient is solved by using the finite element method of a geometry which resembles colonic crypts.
Stochastical time steps are taken via events calculated by Gillespie's algorithm and the resulting model simulates cells wandering along the colonic crypts and releasing when there is enough pressure on a given cell. Biologically, the shedding of the cells are noticed to occur at the surface of the crypts. However, in this model cells release more frequently in the bottom due to the pressure distribution. This could imply that the cells have a more active nature in order to be able to release at the top of the crypts.
ISSN: 1401-5757, TVE-F 17 013 juni Examinator: Martin Sjödin
Ämnesgranskare: Natalia Ferraz Handledare: Stefan Engblom
1 Popul¨ arvetenskaplig sammanfattning
F¨or att f¨orst˚a och ˚aterskapa beteenden av biologiska processer och f¨or¨andringar i m¨anniskokroppen, har det under de senaste ˚artiondena lagts mycket fokus p˚a matematisk modellering. Matematisk modellering anv¨ands idag frekvent inom naturvetenskapliga samt samh¨allsvetenskapliga sammanhang f¨or att f¨orst˚a olika system och f¨oruts¨aga beteenden.
M˚alet med detta projekt ¨ar att ˚aterskapa cellers spridning i den m¨anskliga tarmen. Miljontals sm˚a fack, kallade krypter, utg¨or tarmens geometri, med stamceller som st¨andigt f¨or¨okar sig. Tarmcancer tros ha sitt ursprung fr˚an mutationer inom dessa krypter. Spridningen av celler samverkar mellan krypterna och celler som n˚ar ytan kommer p˚a sl¨appa fr˚an tarmv¨aggen n¨ar cellerna stigit upp ur krypterna.
Med hj¨alp av en modell i programmet MATLAB om cellers spridning kommer interaktionen mellan olika krypter simuleras. Med matematiska metoder l¨oses en tryckf¨ordelning hos cellpopulationer och
slumpm¨assiga h¨andelser driver simulationen. H¨andelserna leder till spridning av de f¨ornyande cellerna i krypterna och om trycket blir nog stort sl¨apper cellerna fr˚an tarmv¨aggen.
Den framtagna modellen klarar av att simulera denna passiva process, men p˚a grund av att trycket p˚a cellerna nere i krypterna ¨ar st¨orre sl¨apper st¨orsta delen celler i botten. Detta kan antyda p˚a att cellspridningen sker via en aktiv mekanism.
Contents
1 Popul¨arvetenskaplig sammanfattning 2
2 Introduction 4
3 Theory 4
3.1 Biological theory . . . . 4
3.2 Mathematical theory . . . . 5
3.2.1 Finite element method . . . . 5
3.2.2 Gillespie’s algorithm . . . . 6
4 Model 6 4.1 Basic model . . . . 6
4.2 Creating geometry . . . . 6
4.3 Initiating the simulation . . . . 7
4.4 Degrees of freedom and pressure gradient . . . . 8
4.5 Probability intesities . . . . 9
4.6 Events . . . . 10
5 Results 10 5.1 Function of model . . . . 10
5.2 Expanded geometry . . . . 11
6 Discussion 12 6.1 General behaviour . . . . 12
6.2 Passive or active process . . . . 13
6.3 Difficulties . . . . 13
7 Conclusion 14
8 Appendix 16
2 Introduction
Mathematical modeling is a powerful tool to investigate hypothesis of di↵erent observed biological, physical or chemical behaviours. One example is to simulate how cells proliferate in the intestine. The process which the model will try to simulate is a close up at the cell proliferation of the outer layer of the intestine.
This layer is built up by creases called colonic crypts and the simulation will investigate the interaction between them. The process is driven by stem cells residing in the bottom of the crypt. Proliferation pushes the cells along the wall and interacts with cells from other crypts due to overpopulation of cells they will shed into the intestine.
Mathematically, this is modelled as each crypt having a single stem cell in the bottom. These cells are the only ones allowed to proliferate. As the stem cell divides, a single slot gains one cell, to then contain two cells and therefore a higher pressure originates from that place. The Laplace-equation is solved in order to gain a pressure gradient of the geometry, which drives to proliferation. Di↵erent probabilities depending on the pressure will force the cells o↵ the boundary. If a cell loses contact with the boundary it will be excluded from the model from that time forward. MATLAB discretisizes the geometry and to solve the pressure gradient, the finite element method is used. The events and the duration of the time step are stochastically calculated with Gillespie’s algorithm. Consequently, the model is event driven.
The aim of the project is to be able to simulate the behavior of cell shedding in the intestine, recreating biological behavior. This is done to expand knowledge of how cells proliferate under specif circumstances.
In order to grasp the di↵erence in behaviour by changing the reaction constants and geometries, some schematic graphic and animations will be included in the report.
3 Theory
3.1 Biological theory
The human intestine contain an outer layer of colonic crypts that is renewed every 2-3 days due to
long-living stem cells that proliferate at the bottom of the crypts.[1] In the large intestine, around 2·107 of tube formed crypts form a creased geometry, see Figure 1.. This is one of the most rapidly renewing mammalian tissues with a complete turnover of 3-5 days. The overpopulation of cells down the crypts builds pressure and the di↵erentiated cells ascend up the crypt wall and as the cells reach the surface they undergo apoptosis, i.e natural cell death, and are shed into the intestinal lumen. [2]
Due to the cell growth in the base of the crypt, the migration is thought to be a mitotic activity. The process is described as a passive pressure-driven process that force the cells up the wall and ultimately release into the lumen. However, an active cell movement with an active base membrane flow and cell shedding is also a possibility. [1]
Figure 1: Fluorescence microscopy image of colonic crypts inside the intestine. The dark parts represent the crypts.[3]
3.2 Mathematical theory
In the simulation two main computational method are used. The finite element method is used to solve the pressure gradient of the geometry and Gillespie’s algorithm is used to advance the proliferation.
3.2.1 Finite element method
Finite element method, called FEM, is a numerical method often used for big simulations of mechanical engineering, biomechanics, aeronautical etc. The simulations often solves problems of heat transfer, pressure gradients or electromagnetic potential. One of the biggest strengths of the FEM-method is its potential to work with complex geometries.
The first step of the method involves discretizing the space to find a solution to a partial di↵erential equation. This is often done by triangulation, especially for 2-dimensional problems. It is also possible to use polygons or squares. The discretization results in a mesh consisting of nodes. From this, a basis function is chosen. Commonly piece-wise smooth and linear. To each node there is a base function, i, assigned. The base functions are continuous throughout the space but are only nonzero in the region surrounding the node. An approximation of the solution, Uh, is solved by computing a weight function uj
for each node. Such that the approximation can be solved with a linear system of equations, in the form :
Uh= Xn j=1
uj j
where n is the number of nodes.[4]
3.2.2 Gillespie’s algorithm
The Gillespie algorithm is an algorithm within probability theory that generates a statistically correct solution of a stochastic equation. A short summary of the steps, somewhat modified to the process used, are presented below:
Step 1 Initialize cell population and define reaction constants for di↵erent events.
Step 2 Generate random number to determine which event that will occur, as well as the time interval to the next event. The probability of an event to occur is proportional to the number of possible events . Step 3 Increase the simulation time with the time step randomly generated in Step 2. Depending on what
event that occurred, update cell population.
Step 4 Unless simulation time is expired, go back to Step 2. [5]
4 Model
The model that will simulate the proliferation of cells is driven by Gillespie’s algorhitm. Consequently, the simulation is event driven, such that the time steps are depending on the number of possible events and not by a given time interval. Each step solves a pressure gradient with the Finite element method, through discretization with triangulation, and a system of linear equations with matrices derived from the method.
Following this short introduction is a closer look at the di↵erent parts of the model.
4.1 Basic model
An already existing model which simulates cells proliferating out through a simple geometry is used to build the model. It has both the Gillespie algorithm implemented and the FEM solution for the pressure gradient.
The basic model has a working visualization using triangulation of the points from the mesh. From the triangulation the model uses a method to create voxels. These voxels are in one of three di↵erent states.
Empty, singly occupied or doubly occupied and the states are represented with di↵erent colours.
4.2 Creating geometry
The first step of the model is creating a geometry which resembles the colonic crypts. It is done with MATLAB’s PDE toolbox. The crypts are simplified to triangles with ratios of depth, width and length similar to actual colonic crypts. [6]. The line between the crypts represents the intestine wall. This is also referred to as the top of the crypts. The second step of the finite element method is triangulation which is also done using PDE toolbox, see Figure 2. This results in a mesh of the geometry, where the parameters that represents the triangulation, in form of vertices, points and edges, are exported from the toolbox and
used in the simulation. Each node of the triangulation represents a slot where cells can exist. Figure 2 shows a geometry of two colonic crypts with the triangulation from PDE toolbox.
Figure 2: Geometry resembling the colonic cryspts cre- ated by Matlabs PDE tool box. The surface of the crypts is represented by y-axis=0, and the bottom of the crypt is everything below. The area above the y-axis represents the lumen were cells can shed when released from the boundary.
4.3 Initiating the simulation
The FEM step of creating a linear system of equations is done by exporting parameters from the mesh of PDE toolbox. From the parameters a sti↵ness matrix and a mass matrix is created. These matrices are used to solve the pressure gradient of the system.
The initial population is created before the simulation begins, which is the first step of the Gillespie algorithm. The nodes in the bottom of each crypt are the only ones which are allowed to proliferate, i.e.
stem cells.
A quotient between the length and edges is also needed to be able to choose events correctly, called gradquotient. The gradquotient is an integration of the edge lengths over the distance from each node in a given direction, shown in figure 3. This is computed using variables from the voxels created in the basic model 4.1.
Figure 3: Edge length and the distance from each node that are used to calculate the gradquotient. The red line is the distance, whereas the highlighted edge is the edge length.
4.4 Degrees of freedom and pressure gradient
To classify and keep track of the state of each voxel, di↵erent ”Degree of freedom” are created to represent the event that may occur. All active cells as well as the empty voxels, which are located next to active cells, form the base of the model. The singly occupied voxels represents the boundary movement and how the population may spread. The overly occupied voxels represents how the pressure builds in the geometry, forcing the cells to spread. Finally, the layer above of the boundary represents the cells that release from the geometry due to pressure. Cells that are released from the geometry are considered dead and removed from the geometry in the next visualization frame. The states are enumerated in the discretized geometry to determine the the events. See section 4.6. For a visualization of the states, see Figure 4.
Figure 4: Visualization of cell proliferation of two crypts , that explains the di↵erent states each voxel can represent.
The pressure Laplacian is created from all active cells, however the model is driven by pressure from overly occupied cells. When a cell is moved due to pressure and the movement changes the geometry, the sti↵ness matrix is updated to get the correct pressure Laplacian. This is used to update the pressure gradient.
When proliferation or movement between already occupied cells occurs, an update of the sti↵ness matrix is not required but the pressure gradient is updated.
4.5 Probability intesities
The second step of Gillespie’s Algorithm, i.e. chose a random event, is taken via computing the probability intensity of either a boundary movement, movement of a cell from a doubly occupied voxel or a cell division event. The movement intensities are proportional to the number of cells in the di↵erent states, but also through the reaction constants previously mentioned. The reaction constants scale the intensities of di↵erent events, allowing some events to happen more frequently, they are also used to modify the simulation to a desired behaviour. The proliferation rate is proportional to the number of voxels in the simulation.
The event and the time step for the next event are both calculated stochastically, according to Gillespie’s algorithm. Both calculations are scaled to the number of possible events and the values of the probabilities, such that if many events have a high probability to happen the duration of the time step decreases.
4.6 Events
There are three possible events. Singly occupied voxels, representing the boundary movement, can move to an empty voxel. The doubly occupied voxels are able to move one cell to either a empty voxel or to a nearby singly occupied voxels. The third possible event is cell division, where a singly occupied voxel becomes doubly occupied.
Using the ”Degrees of Freedom”, presented in 4.4, the movements are done by first finding possible voxels to move from. Given a possible event of a voxel, there can be several di↵erent neighbors to migrate to. Due to the pressure gradient and reaction constants, there are di↵erent probabilities to where it actually moves, and ultimately this event is determined stochastically.
For visualization purposes, the cell population’s state are saved with predetermined time steps depending on the duration of the simulation. Thus, numerous events occurs between every visualization step.
5 Results
5.1 Function of model
By creating geometries that fulfill the biological requirements and by triangulating with a spacing that is not too refined, the resulting meshes according to figure 2 have good conditions for the objective. Through adjusting the reaction constants, the model’s cells are initially able to ”wander” up along the crypt’s walls.
After a layer of single cells have filled the voxels of the crypt, doubly occupied voxels starts appearing on the walls of the crypts. Eventually, the surface gets overpopulated as well, forcing the cells to release from the boundary. A ratio of 25% of cells that were shed at the surface of the geometry were achieved.
Figure 5: Discrete frames of the simulation, showing the behaviour of the model.
5.2 Expanded geometry
An expanded geometry with four proliferating crypts was simulated. Using the same simulation time, the ratio of cells that were shed on the surface were 25% and due to larger probability intensities the time steps were shorter and thus more events occurred.
Figure 6: Simulation image using a geometry with four crypts.
6 Discussion
6.1 General behaviour
The model used in this simulation has a relatively simple approach. Despite this, by scaling the reaction constants correctly, the desired behavior was somewhat achieved. The desired behaviour was to make the cells release at the surface rather than at the bottom and his was achieved when the movement between voxels on the boundary were more likely to occur rather than releasing from the boundary. If the constants are inadequate, all cells release in the bottom due to the pressure always being greater than at the surface of the geometry. Despite this, with appropriate reaction constants the model achieved 25% ratio of cells that were shed at the surface. To further improve this, the reaction constant that gives a cell shedding event could be altered, making it more likely to release at the surface rather than releasing in the crypt. This method would however imply that the cells ”know” their position and therefore suggest an active process.
Comparing the di↵erent geometries, two and four crypts, some logical behavior is observed. Both
geometries builds larger pressure at the center, resulting in more overly occupied voxels with more released cells in the center region. Due to the edges of the geometries being open and not a↵ected by pressure from their neighbouring crypt, the important simulation behavior is limited to the center. Another noticeable behavior is the case of more than two crypts, the whole geometry acts as a pressure gradient with a higher pressure at the center of the population, as seen in figure 6. Red voxels represents an active pressure and a higher density of red voxels indicates that there is an overall greater pressure in that area. As the density of red voxels decreases towards the edges the overall pressure diminishes as well. All this could be translated as a macro pressure gradient, with its greatest pressure at the center.
Because of the model being event driven another unwanted behaviour appears more clearly as the
geometry increases in size. Several events cannot happen at the same time. This interferes with a constant proliferation of stem cells or cells moving far from each other.
6.2 Passive or active process
The cell density in the crypt exceeds the density of the surface of the geometry. This leads to a higher pressure in the crypts, resulting in more cells that release in the bottom. According to the biological theory, the release will mainly occur after the cells ascended to the surface. If the process is passive, the rates of the proliferation is presumably di↵erent than our model, forcing the cells upward the boundary rather than releasing at the bottom. Possible solutions to improve the model would be to track the age of the cells and making them more likely to release after a certain age. Another possible improvement is to alter the reaction constants such that they depend on a acidity gradient. The intestine would vary in acidity depending on the depth, making the cells more likely to release at the surface.
Due to the fact that there always will be a higher pressure in the bottom of the crypt, the likelihood for cells to release at the bottom will always be greater than for them to shed at the top of the geometry.
Because of this inherent pressure di↵erence it would be hard to get the correct behaviour and it is reasonable that it could be an active process.
6.3 Difficulties
One problem throughout this project was the visualization aspect. The visualization from the basic model only works on a convex geometry as the method creates open voxels with extended lines to infinity at the edge of the geometry. On a non convex geometry, these lines intersect and creates misshaped voxels, see Figure 7. As the geometry was non convex this problem was encountered. This was solved by scaling down the visualization to a more simple approach. The voxels are still a part of the code via calculating
gradquotients, however in the visualization only the nodes are highlighted to represent the state each voxel is in.
Since the gradquotient is calculated with values from the voxels, some edge lengths were incorrect. The method creates abnormal voxels with incorrect edge lengths, which leads to some inaccurate probability intensities. Also as mentioned before, some lines are of infinity length and when there is no pressure di↵erence, it leads to Not-a-Number values in the probability intensities.
Figure 7: Convex and non convex voxel geometries displaying the problem of misshaped voxels and incorrect edge lengths.
7 Conclusion
A few behaviours are apparent in the model, mainly that the pressure is at its greatest in the bottom of the crypts. Ergo more cells will always shed in the crypts instead of the surface. This is an inherent problem of the model. It could be improved by implementing reaction constants that changes with the age of the cells, a pH gradient over the crypts or improving the geometry to more accurately resemble the crypts. However, regardless of the improvements the inherent nature of the pressure gradient implies that the biological process might be an active mechanism. It is also apparent that bigger geometries with more crypts makes the model crackle, in form of a macro pressure gradient and events not being able to happen independently of each other. Therefore the interaction of the colonic crypts simulated with this model is most accurately depicted with only two crypts.
References
[1] Matthew D. Johnston, Carina M. Edwards, Walter F. Bodmer, Philip K. Maini, and S. Jonathan Chapman. Examples of mathematical modeling – tales from the crypt. Cell Cycle 2007, 6:2106-2112.
[2] Gary R. Mirams, Alexander G. Fletcher, Philip K. Maini, Helen M. Byrne. A theoretical investigation of the e↵ect of proliferation and adhesion on monoclonal conversion in the colonic crypt. Journal of Theoretical Biology 2012, 312:143-156.
[3] Uppsala University Picture acquired from Mia Phillipsson, Medical Cell Biology. May 2017.
[4] Finite Element Method, https : //en.wikipedia.org/wiki/F inite element method. Last update: May 10, 2017.
[5] Gillespie algorithm, https : //en.wikipedia.org/wiki/Gillespie algorithm. Last update: May 2, 2017.
[6] Matthew D. Johnston, Carina M. Edwards, Walter F. Bodmer, Philip K. Maini, and S. Jonathan Chapman. Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer. PNAS 2007; 104:4008-4013.
8 Appendix
1 %P r o l i f i r a t i o n o f two stem c e l l s i n s e p a r a t e c r y p t s .
2 % Two c r y p t s
3
4 %Martin Edin , N i l s E r l a n s o n
5
6 %r e p e a t a b l e r e s u l t s
7 rng ( 2 2 2 ) ;
8
9 %C r e a t i n g mesh , r e c t a n g u l a r geometry with two c r y p t s .
10 %Loads PDEmesh with output v a r i a b l e s p , e and t .
11 l o a d(’ Geomesh5 . mat ’) ;
12
13 %T r i a n g u l a t e , used t o c a l c u l a t e a g r a d q u o t i e n t .
14 DT = d e l a u n a y T r i a n g u l a t i o n ( p ’ , e ( 1 : 2 , : ) ’ ) ;
15 [ V, R ] = voronoiDiagram (DT) ;
16
17 % s w i t c h format somewhat
18 patchmax = max( c e l l f u n (’ p r o d o f s i z e ’, R ) ) ;
19 R = nan (s i z e( R , 1 ) , patchmax ) ;
20 f o r i = 1 :s i z e(R, 1 ) , R( i , 1 : s i z e( R { i } , 2 ) ) = R { i } ; end
21 % ( a l l o w s f o r v i s u a l i z a t i o n by t h e patch f u n c t i o n )
22
23 % each node i n t h e mesh i s a ( midpoint o f a ) v o x e l
24 Nvoxels = s i z e( p , 2 ) ;
25
26 % a s s e m b l e minus t h e L a p l a c i a n on t h i s g r i d ( i g n o r i n g BCs ) a s w e l l a s
27 % t h e Mass matrix
28 [ L ,M] = assema ( p , t , 1 , 1 , 0 ) ;
29
30 % t h e ( lumped ) mass matrix g i v e s t h e e l e m e n t volume
31 dM = f u l l(sum(M, 2 ) ) ;
32
33 % e x p l i c i t l y i n v e r t t h e lumped mass matrix and f i l t e r t h e d i f f u s i o n matrix
34 [ i , j , s ] = f i n d(L) ;
35 s = s . /dM( i ) ;
36 keep = f i n d( i ˜= j ) ;
37 i = i ( keep ) ; j = j ( keep ) ; s = s ( keep ) ;
38
39 % r e b u i l d L , e n s u r i n g t h a t t h e d i a g o n a l e q u a l s minus t h e sum o f t h e
40 % o f f d i a g o n a l e l e m e n t s
41 L = f s p a r s e ( i , j , s , [ Nvoxels Nvoxels ] ) ;
42 L = L+f s p a r s e ( 1 : Nvoxels , 1 : Nvoxels , f u l l(sum(L , 2 ) ) ’ ) ;
43
44 % s t a t i c n e i g h b o r matrix
45 N = f s p a r s e ( i , j , 1 ,s i z e(L) ) ;
46 %Neighbor r e p r e s e n t a t i o n f o r each node
47 n e i g h = f u l l(sum(N, 2 ) ) ;
48
49 % d i s t a n c e between two v o x e l m i d p o i n t s : D1( k ) = norm (P ( : , i ( k ) ) ,P ( : , j ( k ) ) )
50 D1 = s q r t(sum( ( p ( : , i ) p ( : , j ) ) . ˆ 2 ) ) ’ ;
51
52 % edge l e n g t h : L1 = edge l e n g t h o f edge between v o x e l i and j
53 L1 = z e r o s(s i z e(D1) ) ;
54 f o r k = 1 : numel ( i )
55 n = f s e t o p (’ i n t e r s e c t ’,R( i ( k ) , : ) ,R( j ( k ) , : ) ) ;
56 n (i s n a n( n ) ) = [ ] ;
57 L1 ( k ) = norm(d i f f(V( n , : ) ) ) ;
58 end
59
60 % G r a d q u o t i e n t
61 g r a d q u o t i e n t = f s p a r s e ( i , j , L1 . / D1 ,s i z e(L) ) ;
62
63 %L o c a t i n g t h e two stem c e l l s i n bottom o f c r y p t s
64 [ ˜ , i 1 ] = min( ( p ( 1 , : ) +0.75) .ˆ2+( p ( 2 , : ) +1) . ˆ 2 ) ;
65 [ ˜ , i 2 ] = min( ( p ( 1 , : ) 0.75) .ˆ2+( p ( 2 , : ) +1) . ˆ 2 ) ;
66
67 %S p e c i f i e d c e l l s t h a t p r o l i f i r a t e
68 pdof = [ i 1 i 2 ] ’ ;
69 70 71
72 %l o c a t i n g boundary
73
74 bound = boundary ( p ’ , 1 ) ;
75 keep = f i n d( p ( 2 , bound ) < 0 . 0 2 ) ;
76 bound = bound ( keep ) ;
77 78
79 %L o c a t i n g t h e l a y e r ” j u s t above ” t h e i n t e s t i n e boundary , c a l l e d death DOF(
ddof )
80 %S p e c i f i e d l a y e r o f DOFs where U=0. C e l l s t h a t move t o t h i s DOF r e l e a s e s . .
81 %. . from t h e i n t e s t i n e w a l l and i s c o n s i d e r e d dead .
82 U = f s p a r s e ( bound , 1 , 1 , [ Nvoxels 1 ] ) ;
83 84
85 ddof = f i n d(N⇤(U ˜= 0) > 0 & U == 0) ;
86
87 % C r e a t i n g a boundary matrix t o a c q u i r e d i f f e r e n t r e a c t i o n c o n s t a n t s f o r
88 % boundary movement
89 VU = (U > 0 ) ;
90 91
92 %C r e a t i n g i n i t i a l p o p u l a t i o n c o n s i s t i n g o f two stem c e l l s i n t h e bottom o f
93 %t h e c r y p t s
94
95 U = f s p a r s e ( pdof , 1 , 1 , [ Nvoxels 1 ] ) ;
96
97 % d i f f u s i v e p r e s s u r e r a t e s
98 Drate1 = 0 . 0 0 0 1 ; % from boundary t o non boundary
99 Drate2 = 1 0 0 0 ; % boundary > boundary
100 Drate3 = 0 . 0 1 ; % i n t o a l r e a d y o c c u p i e d v o x e l
101 D r a t e = [ Drate1 Drate2 ; NaN Drate3 ] ;
102
103 %p r o l i f e r a t i o n r a t e
104 lam=2/ Nvoxels ;
105
106 % s i m u l a t i o n i n t e r v a l
107 Tend = 3 0 0 0 0 ;
108 % s o l u t i o n r e c o r d e d a t t h i s g r i d :
109 t s p a n = l i n s p a c e( 0 , Tend , 2 0 1 ) ;
110
111 %S p a r s e s o l u t i o n v e c o r D t o r e c o r d c e l l s t h a t r e l e a s e from t h e boundary t o
112 %ddof
113 D=z e r o s( Nvoxels , 1 ) ;
114 D=f s p a r s e (D) ;
115
116 % r e p r e s e n t a t i o n o f s o l u t i o n : c e l l v e c t o r o f s p a r s e m a t r i c e s
117 %Usave w i l l r e p r e s e n t a c t i v e c e l l s , Dsave w i l l r e p r e s e n t dead c e l l s .
118 Usave = c e l l ( 1 , numel ( t s p a n ) ) ;
119 Dsave = c e l l ( 1 , numel ( t s p a n ) ) ;
120
121 i f ˜e x i s t(’ r e p o r t p r o g r e s s ’,’ var ’)
122 r e p o r t p r o g r e s s = t r u e ;
123 end
124 i f r e p o r t p r o g r e s s
125 r e p o r t ( tspan , U, ’ i n i t ’) ;
126 e l s e
127 r e p o r t ( tspan , U, ’ none ’) ;
128 end
129
130 t t = t s p a n ( 1 ) ;
131 Dsave{1} = D;
132 Usave{1} = U;
133 i = 1 ;
134 % l o g i c f o r r e u s e o f LU f a c t o r i z a t i o n s
135 updLU = t r u e ;
136 La = s t r u c t (’X ’, 0 , ’L ’, 0 ,’U ’ , 0 ,’ p ’ , 0 ,’ q ’, 0 , ’R ’, 0 ) ;
137
138 %c o u n t e r t o check e v e n t s w i t h i n w h i l e l o o p and number o f dead c e l l s
139 c o u n t e r =0;
140 p c o u n t e r =0;
141 t w o c o u n t e r =0;
142 o n e c o u n t e r =0;
143 d e a t h c o u n t e r U P =0;
144 death counter DOWN =0;
145
146 w h i l e t t <= t s p a n (end)
147
148 % c l a s s i f y t h e Degrees Of Freedoms (DOFs f o r s h o r t )
149
150 % a c t i v e DOFs : o c c u p i e d v o x e l s
151 a d o f = f i n d(U) ;
152 % boundary DOFs which may move : c o n t a i n i n g one c e l l p e r v o x e l and
153 % with an empty v o x e l b e s i d e s i t
154 bdof m = f i n d(N⇤(U > 0) < neigh & U == 1) ;
155 % s o u r c e DOFs c o n t a i n i n g more than one c e l l p e r v o x e l ( a c t u a l l y 2 )
156 s d o f = f i n d(U > 1 ) ;
157
158 % s o u r c e DOFs which may move : with a v o x e l c o n t a i n i n g l e s s number o f
159 % c e l l s next t o i t ( a c t u a l l y 1 o r 0 )
160 sdof m = f i n d(N⇤(U > 1) < neigh & U > 1) ;
161
162 % i n j e c t i o n DOFs : t h e l a y e r o f v o x e l s ” j u s t o u t s i d e ” a d o f where we may
163 % i n j e c t a BC ( p r e s s u r e 0 o f t h e f r e e matrix )
164 I d o f = (N⇤(U ˜= 0) > 0 & U == 0) ;
165 i d o f = f i n d( I d o f ) ;
166
167 % ” A l l DOFs” = a d o f + i d o f , l i k e t h e ” h u l l o f a d o f ”
168 Adof = [ a d o f ; i d o f ] ;
169
170 % The above w i l l be enumerated w i t h i n U, a Nvoxels by Nvoxels s p a r s e
171 % matrix . Determine a l s o a l o c a l enumeration , eg . [ 1 2 3
172 % . . . numel ( Adof ) ] .
173 Adof = ( 1 : numel ( Adof ) ) ’ ;
174 [ ˜ , i x ] = f s e t o p (’ ismember ’, bdof m ’ , Adof ’ ) ;
175 bdof m = Adof ( i x ) ;
176 [ ˜ , i x ] = f s e t o p (’ ismember ’, s d o f ’ , Adof ’ ) ;
177 s d o f = Adof ( i x ) ;
178 [ ˜ , i x ] = f s e t o p (’ ismember ’, sdof m ’ , Adof ’ ) ;
179 s d o f m = Adof ( i x ) ;
180 [ ˜ , i x ] = f s e t o p (’ ismember ’, i d o f ’ , Adof ’ ) ;
181 i d o f = Adof ( i x ) ;
182 [ ˜ , i x ] = f s e t o p (’ ismember ’, pdof ’ , Adof ’ ) ;
183 p d o f = Adof ( i x ) ;
184 185
186 i f updLU
187 % p r e s s u r e L a p l a c i a n f a c t o r i z a t i o n i s r e u s e d
188 La .X = L( Adof , Adof ) ;
189
190 % s e l e c t i n g a l l i n j e c t i o n DOFs
191 Lai = f s p a r s e ( i d o f , i d o f , 1 ,s i z e( La .X) ) ;
192 % remove e q s ( rows ) f o r i n j e c t i o n DOFs, r e p l a c e with d i r e c t i n j e c t i o n
193 La .X = La . X Lai⇤La .X+Lai ;
194
195 % f a c t o r i z e
196 [ La . L , La . U, La . p , La . q , La .R] = l u( La . X, ’ v e c t o r ’) ;
197 updLU = f a l s e ; % assume we can r e u s e
198 end
199
200 % RHS s o u r c e term p r o p o r t i o n a l t o t h e over occupancy ( t h e r e s t o f t h e
201 % z e r o s e n s u r e t h e homogeneous D i r i c h l e t BC i s s a t i s f i e d a t i d o f )
202
203 Pr = f u l l( f s p a r s e ( s d o f , 1 , 1 . /dM( s d o f ) , [s i z e( La . X, 1 ) 1 ] ) ) ; % RHS f i r s t . . .
204 Pr ( La . q ) = La .U\( La . L\( La .R( : , La . p ) \Pr ) ) ; % . . then t h e s o l u t i o n
205
206 % compute i n t e n s i t i e s o f p o s s i b l e e v e n t s
207
208 % moving boundary DOFs : each empty n e i g h b o u r v o x e l i s a s s o c i a t e d with
209 % a f l o w r a t e p r o p o r t i o n a l t o t h e p r e s s u r e g r a d i e n t i n t h a t
210 % d i r e c t i o n i n t e g r a t e d o v e r t h e c o r r e s p o n d i n g edge
211
212 [ i i , j j , gq ] = f i n d( g r a d q u o t i e n t ( bdof m , Adof ) ) ; % n e i g h b o u r s . . .
213 keep = f i n d(U( Adof ( j j ) ) == 0 ) ; % . . . t o move t o
214 i i = r e s h a p e( i i ( keep ) , [ ] , 1 ) ;
215 j j = r e s h a p e( j j ( keep ) , [ ] , 1 ) ;
216 gq = r e s h a p e( gq ( keep ) , [ ] , 1 ) ;
217 grad = f s p a r s e ( i i , 1 , gq .⇤max( ( Pr ( bdof m ( i i ) ) Pr ( j j ) ) , 0 ) . . .
218 .⇤ Drate (2⇤VU( Adof ( j j ) ) +1) , numel ( bdof m ) ) ;
219 moveb = f u l l( grad ) ;
220
221 % c h e c k i n g f o r c a s e s where moveb g i v e s i n f⇤0 and removing that event
222 i f any(i s n a n( moveb ) )
223 check = i s n a n( moveb ) ;
224 moveb ( check ) = 0 ;
225 end
226 227 228
229 % C e l l s i n a doubly o c c u p i e d v o x e l may move by t h e same p h y s i c s
230 [ i i , j j , gq ] = f i n d( g r a d q u o t i e n t ( sdof m , Adof ) ) ; % n e i g h b o u r s . . .
231 keep = f i n d(U( Adof ( j j ) ) < 2 ) ; % . . . t o move t o
232 i i = r e s h a p e( i i ( keep ) , [ ] , 1 ) ;
233 j j = r e s h a p e( j j ( keep ) , [ ] , 1 ) ;
234 gq = r e s h a p e( gq ( keep ) , [ ] , 1 ) ;
235 % remove any p o s s i b l y r e m a i n i n g n e g a t i v e r a t e s
236 grad = f s p a r s e ( i i , 1 , gq .⇤max( Pr ( s d o f m ( i i ) ) Pr ( j j ) , 0 ) .⇤ . . .
237 D r a t e ( 2⇤VU( Adof ( j j ) )+U( Adof ( j j ) ) +1) , numel ( sdof m ) ) ;
238 moves = f u l l( grad ) ;
239
240 % c h e c k i n g f o r c a s e s where moves g i v e s i n f⇤0 and removing that event
241 i f any(i s n a n( moves ) )
242 check = i s n a n( moves ) ;
243 moves ( check ) = 0 ;
244 end
245 246
247 % P r o l i f e r a t i o n i n t e s i t y
248 i f any(U( pdof )>=2)
249 p r o l i f =0;
250 e l s e
251 p r o l i f = lam⇤U( pdof ) ;
252 end
253
254 %I n t e n s i t i e s o f p o s s i b l e e v e n t s : Boundary , s o u r c e and p r o l i f i r a t i o n
255 i n t e n s = [ moveb ; moves ; p r o l i f ] ;
256
257 % w a i t i n g time u n t i l next e v e n t and t h e e v e n t i t s e l f
258 lambda = sum( i n t e n s ) ;
259 dt = r e a l l o g (rand) / lambda ;
260 rnd = rand⇤lambda ;
261 cum = i n t e n s ( 1 ) ;
262 i x = 1 ;
263 w h i l e rnd > cum
264 i x = i x +1;
265 cum = cum+i n t e n s ( i x ) ;
266 end
267 % ( now i x p o i n t s t o t h e i n t e n s i t y which f i r e d f i r s t )
268
269 % Record v a l u e s f o r v i s u a l i z a t i o n , Usave f o r a c t i v e c e l l s and Dsave f o r
270 % dead c e l l s .
271 i f t s p a n ( i +1) < t t+dt
272 i e n d = i+f i n d( t s p a n ( i +1:end) < t t+dt , 1 ,’ l a s t ’) ;
273 Usave ( i +1: i e n d ) = {U} ;
274 Dsave ( i +1: i e n d ) = {D} ;
275 i = i e n d ;
276 i f i e n d == numel ( Usave ) , break; end
277
278 end
279
280 i f i x <= numel ( moveb )
281 % movement o f a boundary ( s i n g l y o c c u p i e d ) v o x e l
282 i x = bdof m ( i x ) ;
283 i x = Adof ( i x ) ;
284
285 % s e l e c t n e i g h b o r t o move t o
286 j x = f i n d(N( ix , Adof ) ) ;
287 j x = j x (U( Adof ( j x ) ) == 0 ) ; % o n l y a l l o w moves t o l e s s p o p u l a t e d v o x e l s
288 % t h e p r e s s u r e g r a d i e n t d r i v e s t h e movement
289 r a t e s = f u l l( g r a d q u o t i e n t ( ix , Adof ( j x ) ) ) ’ .⇤max( Pr ( i x ) Pr ( j x ) , 0 ) .⇤ . . .
290 D r a t e ( 2⇤VU( Adof ( j x ) ) +1) ;
291 m = f i n d(cumsum( r a t e s ) > rand⇤sum( r a t e s ) , 1 , ’ f i r s t ’) ;
292 n = Adof ( j x (m) ) ;
293
294 %e x e c u t e e v e n t : move from i x t o n
295 %c h e c k i n g i f a c e l l moves t o ddof , then i t r e l e a s e s from boundary
296 i f any( ismember ( n , ddof ) )
297
298 U( i x ) = U( i x ) 1; %From a c t i v e . .
299 U( n ) =0; %. . t o dead
300 D( n )=D( n ) +1; %C e l l moved t o ddof and r e c o r d e d
301 %Record where t h e c e l l s a r e r e l e a s e d
302 i f p ( 2 , n ) < 0.02
303 death counter DOWN = death counter DOWN +1;
304 e l s e
305 d e a t h c o u n t e r U P = d e a t h c o u n t e r U P +1;
306 end
307 308
309 e l s e %move from i x t o n
310 U( i x ) = U( i x ) 1;
311 U( n ) = U( n ) +1;
312 end;
313 updLU = t r u e ;
314 o n e c o u n t e r=o n e c o u n t e r +1;
315
316 e l s e i f i x <= numel ( moveb )+numel ( moves )
317 % movement o f a c e l l i n a doubly o c c u p i e d v o x e l
318 i x = i x numel ( moveb ) ;
319 i x = s d o f m ( i x ) ;
320 i x = Adof ( i x ) ;
321 % s e l e c t n e i g h b o r t o move t o
322 j x = f i n d(N( ix , Adof ) ) ;
323 j x = j x (U( Adof ( j x ) ) < 2 ) ; % o n l y a l l o w moves t o l e s s p o p u l a t e d v o x e l s
324 % t h e p r e s s u r e g r a d i e n t d r i v e s t h e movement
325 r a t e s = f u l l( g r a d q u o t i e n t ( ix , Adof ( j x ) ) ) ’ .⇤max( Pr ( i x ) Pr ( j x ) , 0 ) .⇤ . . .
326 D r a t e ( 2⇤VU( Adof ( j x ) )+U( Adof ( j x ) ) +1) ;
327 % ( t h i n n i n g o f n e g a t i v e g r a d i e n t s )
328 m = f i n d(cumsum( r a t e s ) > rand⇤sum( r a t e s ) , 1 , ’ f i r s t ’) ;
329 n = Adof ( j x (m) ) ;
330
331 % e x e c u t e e v e n t : move from i x t o n
332 i f U( n ) == 0
333 updLU = t r u e ;
334 end
335
336 %c h e c k i n g i f a c e l l i n a doubly o c c u p i e d v o x e l moves t o ddof
337 %then i t r e l e a s e s from boundary
338 i f any( ismember ( n , ddof ) )
339 U( i x ) = U( i x ) 1; %From a c t i v e . .
340 U( n ) =0; %. . t o dead
341 D( n )=D( n ) +1; %C e l l moved t o ddof and r e c o r d e d
342 %Record where t h e c e l l s a r e r e l e a s e d
343 i f p ( 2 , n ) < 0.02
344 death counter DOWN = death counter DOWN +1;
345 e l s e
346 d e a t h c o u n t e r U P = d e a t h c o u n t e r U P +1;
347 end
348
349 e l s e
350 U( i x ) = U( i x ) 1;
351 U( n ) = U( n ) +1;
352 end;
353 t w o c o u n t e r=t w o c o u n t e r +1;
354 e l s e
355 % o t h e r w i s e t h e stem c e l l s p r o l i f i r a t e
356 i x = i x numel ( moveb ) numel ( moves ) ;
357 i x = p d o f ( i x ) ;
358 i x = Adof ( i x ) ;
359 U( i x ) = U( i x ) +1;
360 p c o u n t e r=p c o u n t e r +1;
361 end
362 t t = t t+dt ;
363 r e p o r t ( t t , U, ’ ’) ;
364 c o u n t e r=c o u n t e r +1;
365 end
366 r e p o r t ( t t , U, ’ done ’) ;
367 % r e t u r n ;
368 % S i m u l a t i o n
369 M = s t r u c t (’ c d a t a ’,{ } ,’ colormap ’ ,{ } ) ;
370 f i g u r e( 1 ) ;
371 s e t(g c f, ’ Renderer ’, ’ o p e n g l ’)
372
373 f o r i = 1 : numel ( Usave )
374 p l o t( p ( 1 , ddof ) , p ( 2 , ddof ) ,’ o ’,’ C o l o r ’ , [ 0 0 0 ] ) ;
375 h o l d on
376
377 a x i s( [ 1 . 5 1 . 5 1 1 ] ) ; a x i s e qual , a x i s o f f
378 i i = f i n d( Usave{ i } == 1) ;
379 p l o t( p ( 1 , i i ) , p ( 2 , i i ) ,’ . ’,’ M a r k e r S i z e ’, 3 0 ,’ C o l o r ’ , [ 0 1 0 ] ) ;
380 j j = f i n d( Usave{ i } == 2) ;
381 p l o t( p ( 1 , j j ) , p ( 2 , j j ) ,’ . ’,’ M a r k e r S i z e ’, 3 0 ,’ C o l o r ’ , [ 1 0 0 ] ) ;
382 p l o t( p ( 1 , pdof ) , p ( 2 , pdof ) ,’ . ’,’ M ar k e r S i z e ’, 3 0 , ’ C o l o r ’ , [ 0 0 1 ] ) ;
383
384 i f i==1
385 kk=f i n d( Dsave{ i } ˜=0) ;
386 p l o t( p ( 1 , kk ) , p ( 2 , kk ) ,’ . ’,’ M a r k e r S i z e ’, 3 0 ,’ C o l o r ’ , [ 0 0 0 ] ) ;
387 e l s e
388 kk=f i n d( Dsave{ i } ˜=Dsave{ i 1}) ;
389 p l o t( p ( 1 , kk ) , p ( 2 , kk ) ,’ . ’,’ M a r k e r S i z e ’, 3 0 , ’ C o l o r ’ , [ 0 0 0 ] ) ;
390 end
391
392 t i t l e(s p r i n t f(’ t = %f ’, t s p a n ( i ) ) ) ;
393 M( i ) = g e t f r a m e(g c f) ;
394 h o l d o f f
395 end
396
397 m o v i e 2 g i f (M,{M( 1 : 2 0 1 ) . cdata } ,’ /home/ n i l s / Documents /MATLAB/Kand . Jobb / s t e n g l i b master / C e l l / a n i m a t i o n s / C r y p t s p r o l i f . g i f ’. . .
398 ,’ d e l a y t i m e ’, 0 , ’ l o o p c o u n t ’, 0 ) ;
