2-dimensional vertex model  simulating epithelial tissue
In[]:=
LaunchKernels[]​​ParallelTable[$KernelID,{i,$KernelCount}]
Out[]=
KernelObject
Name: Local kernel
KernelID: 1
,KernelObject
Name: Local kernel
KernelID: 2
,KernelObject
Name: Local kernel
KernelID: 3
,KernelObject
Name: Local kernel
KernelID: 4
,KernelObject
Name: Local kernel
KernelID: 5
,KernelObject
Name: Local kernel
KernelID: 6

Out[]=
{6,5,4,3,2,1}

geometrical f(x)s

In[]:=
(*getvertexinthelistbyccworder*)​​Clear[getCounterClockwise];​​getCounterClockwise[vertex_,vertices_]:=Block[{pos,v},​​pos=First@@Position[vertices,vertex];​​If[posLength[vertices],pos=1,pos+=1];​​vertices[[pos]]​​];
In[]:=
(*getvertexinthelistbycworder*)​​Clear[getClockwise];​​getClockwise[vertex_,vertices_]:=Block[{ls,pos},​​pos=First@@Position[vertices,vertex];​​If[pos1,pos=Length[vertices],pos-=1];​​vertices[[pos]]​​];
In[]:=
getCounterClockwise[{
x
i
,
y
i
},{{
x
i-1
,
y
i-1
},{
x
i
,
y
i
},{
x
i+1
,
y
i+1
},d,e}]
Out[]=
{
x
1+i
,
y
1+i
}
In[]:=
getClockwise[{
x
i
,
y
i
},{{
x
i-1
,
y
i-1
},{
x
i
,
y
i
},{
x
i+1
,
y
i+1
},d,e}]
Out[]=
{
x
-1+i
,
y
-1+i
}
In[]:=
Clear[areaOfPolygon];​​areaOfPolygon[cells_/;Head[cells]===Association]:=Parallelize[Map[Area@*Polygon,cells]];
In[]:=
Clear[areaPolygon];​​areaPolygon[vertices_]:=Block[{edges},​​edges=Partition[vertices,2,1,1];​​0.5Abs@Total[(#[[1,1]]*#[[2,2]])-(#[[2,1]]*#[[1,2]])&/@edges]​​]
In[]:=
Clear[perimeterOfPolygon];​​perimeterOfPolygon[cells_/;Head[cells]===Association]:=Parallelize[(Perimeter@*Polygon)/@cells];
In[]:=
Clear[perimeterPolygon];​​perimeterPolygon[vertices_]:=Block[{edges},​​edges=Partition[vertices,2,1,1];​​Total[Apply[EuclideanDistance]/@edges]​​]
In[]:=
Clear[centroidPolygon];​​centroidPolygon[vertices_]:=Mean[vertices]
In[]:=
(*counterclockwisepolygonQ*)​​Block{signedarea=0,j,vertlen=5},​​Do[​​j=Mod[i,vertlen]+1;​​signedarea+=(
x
i
y
j
-
x
j
y
i
),​​{i,vertlen}];​​Echo
1
2
(signedarea)​​;
»
1
2
(-
x
2
y
1
+
x
5
y
1
+
x
1
y
2
-
x
3
y
2
+
x
2
y
3
-
x
4
y
3
+
x
3
y
4
-
x
5
y
4
-
x
1
y
5
+
x
4
y
5
)
In[]:=
Clear[polyCounterClockwiseQ];​​polyCounterClockwiseQ[poly_]:=Block[{area=0,j,vertLength=Length[poly]},​​Do[​​j=Mod[i,vertLength]+1;​​area+=poly[[i,1]]*poly[[j,2]];​​area-=poly[[j,1]]*poly[[i,2]],​​{i,vertLength}​​];​​(area/2.)>0​​]
In[]:=
(*Clear[sortCC];​​sortCC[polyinds_,indTopts_,ptsToInds_]:=Block[{cent,poly},​​poly=Lookup[indTopts,polyinds];​​Lookup[ptsToInds,​​DeleteDuplicates@Flatten[MeshPrimitives[ConvexHullMesh[poly],1]/.LineSequence,1]​​]​​];*)
In[]:=
(*sortpointsforaconvexpolygonincounter-clockwisedirection*)
In[]:=
Clear[sortPointsCC];​​sortPointsCC[polyinds_,indTopts_,ptsToInds_]:=Block[{cent,ordering,polyPoints},​​polyPoints=Lookup[indTopts,polyinds];​​cent=Mean[polyPoints];​​ordering=Ordering[ArcTan[#[[1]],#[[2]]]&@(#-cent)&/@polyPoints];​​Lookup[ptsToInds,Part[polyPoints,ordering]]​​]
In[]:=
outeredges[indToPtsAssoc_,localtopology_]:=Block[{k,temp,tcells,assoc},​​Reap[​​Scan[(k=#;temp=localtopology[k];​​tcells=temp[[2]];​​If[tcells≠{},​​MapAt[Sow@Cases[Partition[#,2,1,1],OrderlessPatternSequence[{x:indToPtsAssoc[k],y:_}]{x,y}]&,temp[[1]],{Key[#]}&/@tcells];​​])&,Keys@indToPtsAssoc]​​][[2]]//Flatten[#,2]&​​];

mesh restructuring operations

simulation domain
In[]:=
xLim=yLim={-1,1};​​=Rectangle[{First@xLim,First@yLim},{Last@xLim,Last@yLim}];

local topology

In[]:=
ClearAll[periodicRules,transformRules];​​periodicRules::Information="shift the points outside the simulation domain to inside the domain";​​transformRules::Information="vector that shifts the point outside the simulation domain back inside";​​With[{xlim1=xLim[[1]],xlim2=xLim[[2]],ylim1=yLim[[1]],ylim2=yLim[[2]],dstep=2},​​periodicRules=Dispatch[{​​{x_/;x≥xlim2,y_/;y≤ylim1}SetPrecision[{x-dstep,y+dstep},10],​​{x_/;x≥xlim2,y_/;ylim1<y<ylim2}SetPrecision[{x-dstep,y},10],​​{x_/;xlim1<x<xlim2,y_/;y≤ylim1}SetPrecision[{x,y+dstep},10],​​{x_/;x≤xlim1,y_/;y≤ylim1}SetPrecision[{x+dstep,y+dstep},10],​​{x_/;x≤xlim1,y_/;ylim1<y<ylim2}SetPrecision[{x+dstep,y},10],​​{x_/;x≤xlim1,y_/;y≥ylim2}SetPrecision[{x+dstep,y-dstep},10],​​{x_/;xlim1<x<xlim2,y_/;y≥ylim2}SetPrecision[{x,y-dstep},10],​​{x_/;x≥xlim2,y_/;y≥ylim2}SetPrecision[{x-dstep,y-dstep},10]​​}];​​​​transformRules=Dispatch[{​​{x_/;x≥xlim2,y_/;y≤ylim1}{-dstep,dstep}~SetPrecision~10,​​{x_/;x≥xlim2,y_/;ylim1<y<ylim2}{-dstep,0}~SetPrecision~10,​​{x_/;xlim1<x<xlim2,y_/;y≤ylim1}{0,dstep}~SetPrecision~10,​​{x_/;x≤xlim1,y_/;y≤ylim1}{dstep,dstep}~SetPrecision~10,​​{x_/;x≤xlim1,y_/;ylim1<y<ylim2}{dstep,0}~SetPrecision~10,​​{x_/;x≤xlim1,y_/;y≥ylim2}{dstep,-dstep}~SetPrecision~10,​​{x_/;xlim1<x<xlim2,y_/;y≥ylim2}{0,-dstep}~SetPrecision~10,​​{x_/;x≥xlim2,y_/;y≥ylim2}{-dstep,-dstep}~SetPrecision~10,​​{___Real}{0,0}~SetPrecision~10}];​​];

T1 transition

find edges

T1

T2 transition

find cells for T2

T2

cell division

Difference from Farahdifar's paper: in his article he dilates the cells before dividing them. Here I randomly select a cell and divide it into two. This does not really matter

computing forces

area elasticity

perimeter elasticity

line tension

ΣF

generating mesh

Main( )

miscellaneous