cell_growthΒΆ

../_images/cell_growth.png
#!/usr/bin/env python-mpacts
import mpacts.particles as prt
import mpacts.predicates.cellcycle_predicates as cpreds
import mpacts.predicates.predicates as preds
import mpacts.commands.monitors.progress as pp
from mpacts.core.simulation import Simulation
import mpacts.core.command as cmd
import mpacts.geometrygenerators.quadgeometries as quadgen
import mpacts.contact.matrix.pythontools as cmpt
from mpacts.core.valueproperties import Variable, VariableFunction
from mpacts.commands.misc.stresscalculation import AddStressCalculation
from mpacts.commands.misc.cellcycle import BinaryFissionCommand
from mpacts.commands.misc.cellcycle import Growth_QuCommand
from mpacts.commands.misc.cellcycle import QuiescenceCommand
from mpacts.contact.detectors.multigrid import MultiGridContactDetector
from mpacts.contact.matrix.cmtypes import SumToDiagonalCommand
from mpacts.contact.matrix.cmtypes import ConjugateGradientSolver
from mpacts.contact.matrix.frictionmatrix import FrictionMatrix
from mpacts.contact.models.misc.setmatrix import RemoveResetElements
from mpacts.commands.misc.frictionmatrices import SphereDiffusionMatrixCommand
from mpacts.commands.force.random import ThermalForceCommand
from mpacts.contact.models.collision.celldivision import DividingJKRAreaMatrix
from mpacts.contact.models.collision.jkr.jkr_matrix import JKRMatrix
from mpacts.predicates.shortened_predicates import *
from mpacts.tools.setcontactdatastorage import SetContactDataStorage
from mpacts.core.arrays import create_array
from mpacts.tools.random_seed import set_random_seed
from mpacts.core.units import unit_registry as u
import random

#To ensure that the random numbers are always chosen the same.
#This enables hash checking of the final result (see below)
set_random_seed(1)

#------------------------------------------------------------------------------------------------------------------
#We can use a high timestep of 10 seconds because we solve this as an overdamped system.
#The downside is that one timestep usually takes longer to calculate
mysim = Simulation("simulation", timestep=10*u('s') )
SetContactDataStorage("Vector")
params = mysim('params')

Ecell      = Variable( 'Ecell'        , params, value=1*u("kPa")       )
Ewall      = Variable( "Ewall"        , params, value=1.5*u("kPa")     )
nucell     = Variable( "nucell"       , params, value=0.4              )
nuwall     = Variable( "nuwall"       , params, value=0.4              )
adhcell    = Variable( "adhesion"     , params, value=0.1*u('nN/um')   )
gtime      = Variable( "growth_time"  , params, value=5*u('hr')        )
dtime      = Variable( "division_time", params, value=45*u('min')      )
max_stress = Variable( "max_stress"   , params, value=-5*u('Pa')       )
Rcell      = Variable( "R"            , params, value=6*u('um')        )
Rstd       = Variable( "Rstd"         , params, value=0.2*u('um')      )
gamma      = Variable( "gammacell"    , params, value=2*u('kPa*s/um')  )
gamma_div  = Variable( "gammadiv"     , params, value=500*u('nN*s/um') )
teff       = Variable( "T_eff"        , params, value=310*u('K')       )
visc       = Variable( "visc_eff"     , params, value=0.5*u('Pa*s')    )
cd_ue      = Variable( "cd_ue"        , params, value=10               )
cg_tol     = Variable( "cg_tol"       , params, value=1e-4             )
out_int    = Variable( "out_int"      , params, value=20*u('min')      )
kdiv       = Variable( "k_div"        , params, value=2*u('nN/um')     )
endtime    = Variable( "endtime"      , params, value=3*u('day')       )

R0      = VariableFunction("R0"     , params, function = '0.75*$R$'    )
maxsize = VariableFunction("maxsize", params, function = '1.1*$R$'     )
cd_kd   = VariableFunction("cd_kd"  , params, function = '0.5*$R$'     )

#Make the cells particle container cells:

cells = prt.ParticleContainer('cells',prt.Sphere0, mysim)

#For cells, we need to keep some additional arrays:
create_array("Index","state", cells)
create_array("Scalar","r0", cells)
create_array("Scalar","maxsize", cells)
create_array("Scalar","division_time", cells)
create_array("int","division_partner", cells)
create_array("SymmetricMatrix", "ContactMatrixDiagonal", cells)
create_array("SymmetricMatrix", "chi", cells )
create_array("SymmetricMatrix", "D_sqrt", cells )
create_array("Vector", "div_D_F", cells )

SphereDiffusionMatrixCommand("ComputeFrictionMatrixCmd", mysim
                            , pc = cells, eta = visc, T = teff
                            , gamma = cells['ContactMatrixDiagonal'] )

ThermalForceCommand("ThermalForceCmd", mysim, pc = cells )

#We calculate a (virial) stress-tensor in the contact forces.
#Its eigen values give the principal stress values
AddStressCalculation( mysim, cells, evals = True, evecs = False)

box = prt.ParticleContainer("box",prt.RigidBody.compose( prt.RigidQuad ), mysim)

box("RigidQuads").QuadNormalsCmd( executeOnce=True, x = box('controlPoints')['x'] )

#------------------------------------------------------------------------------------------------------------------
#Cell Cycle commands:

#Define state predicates for growth, division and quiescence
isGrowing = preds.Predicate_StateSaver("isGrowingPredicate", state_idx = 0)
isDividing = preds.Predicate_StateSaver("isDividingPredicate", state_idx = 1)
isQuiescent = preds.Predicate_StateSaver("isQuiescentPredicate", state_idx = 2)

#Define conditions signaling end of growth and end of division. Queiscence never ends ;-)
doneGrowing = cpreds.Predicate_GrowthQuFinished("doneGrowing", predicate_running = isGrowing )

doneDividing = cpreds.Predicate_BinaryFissionFinished( "doneDividing"
                                                     , predicate_running = isDividing
                                                     , child_predicate = isGrowing
                                                     , split_cutoff = 0.95 )

#Condition for high stress:
hasHighStress = preds.Predicate_Stress( "StressPredicate", max_stress = max_stress, component = -1)

#If done growing and high stress, we become quiescent:
quiescenceCondition = AND( doneGrowing, hasHighStress, do = isQuiescent)
#If done growing and not elegible for quiescence, we start dividing:
dividingCondition = AND( doneGrowing, NOT(isQuiescent), do = isDividing)


#Define and init method (based on predicates) for how to initialize cytokinesis. This one is a random direction
div_init_method = cpreds.InitBinaryFission_Random("Init_division", init_offset=1e-3)

#Finally the actual cell cycle commands:

#1) growth command. Start predicate is doneDividing
growthCmd = Growth_QuCommand( "GrowthCommand", mysim, pc=cells
                            , growth_time = gtime, radius_max = maxsize, radius_sd = Rstd
                            , density = 1, max_stress = 1e8, predicate = doneDividing )

#2) Quiescence command: start predicatee is quiescenceCondition
quiescenceCmd = QuiescenceCommand("Quiescence", mysim, pc=cells
                                 , predicate = quiescenceCondition )

#3) Division command: start predicate is dividingCondition
divisionCmd = BinaryFissionCommand( "BinaryFission", mysim, pc=cells
                                  , predicate = dividingCondition
                                  , init_method = div_init_method )
#------------------------------------------------------------------------------------------------------------------
# Contact model - including contact 'division force' pushing the cells aparts
CM_cell_cell = DividingJKRAreaMatrix( pc1=cells, pc2=cells
                                    , E1=Ecell, E2=Ecell
                                    , nu1=nucell, nu2=nucell
                                    , attrConst =adhcell
                                    , division_time = dtime
                                    , division_force = 1#AKA unlimited (1N for cells is huge)
                                    , state_dividing = isDividing.get('state_idx')
                                    , divstiffrep = kdiv
                                    , divstiffatt = kdiv
                                    , gamma_normal = gamma
                                    , gamma_tangential = gamma
                                    , gamma_dividing = gamma_div )

CD_cell_cell = MultiGridContactDetector("CD_cell_cell", mysim
                                       , cmodel=CM_cell_cell
                                       , update_every=cd_ue
                                       , keep_distance = cd_kd )
#Contact models
CM_box_cell = JKRMatrix( pc1=box('RigidQuads'), pc2=cells
                       , E1=Ewall, E2=Ecell
                       , nu1=nuwall, nu2=nucell
                       , attrConst=adhcell
                       , gamma_normal = gamma
                       , gamma_tangential = gamma
                       )

CD_box_cell = MultiGridContactDetector("CD_box_cell", mysim
                                      , cmodel=CM_box_cell
                                      , update_every = cd_ue
                                      , keep_distance = cd_kd )

#------------------------------------------------------------------------------------------------------------------
#make ContactMatrix
CMatrix = FrictionMatrix("ContactMatrix_cells",mysim
                        , pc_list = [ cells ]
                        , cd_list = [ CD_cell_cell ]
                        )
#make command for summing to diagonal
SumDiag_cmd = SumToDiagonalCommand( "SumToDiagonalCmd_cells"
                                  , mysim, A=CMatrix
                                  , Op=cmpt.ContactMatrixOperatorsR3() )

#make ConjugateGradientSolver to update the velocities
ConjugateGradient = cmpt.make_CG_solver( CMatrix, sim=mysim, tolerance=cg_tol, reset_x=False )

#time integration
cells.TimeIntegration_ForwardEuler_Generic_Cmd()

#------------------------------------------------------------------------------------------------------------------
#Add an output writer command for the cells:
cells.VTKWriterCmd(executeInterval = out_int, select_all=True)
box.VTKWriterCmd(executeOnce=True)

#------------------------------------------------------------------------------------------------------------------
#Add a command that prints some information to the screwen at a certain time interval, so we know what we are doing
printer_list = [ pp.DefaultProgressPrinter(mysim,sim_time_unit='h')
               , pp.ManagerSizePrinter(cells)]

printer = pp.PrinterChain(printer_list)

pp.ProgressIndicator("PrintProgress", mysim, printer, print_interval=2.0)

#------------------------------------------------------------------------------------------------------------------

xs = [( 0,-0.99*Rcell.get_value(),0)
      ,(0, 0.99*Rcell.get_value(),0) ]

box_size = [ 6e-5, 6e-5, 6e-5 ]
rs = [random.gauss( Rcell.get_value(), Rstd.get_value() ) for _ in xs]
cells.add_and_set_particles( x = xs, r = rs
                           , state = 0, r0 = R0.get_value()
                           , maxsize = maxsize.get_value(), division_time = -1
                           , division_partner = -1 )

(points, quads) = quadgen.create_box( box_size[0], box_size[1], box_size[2] )

abox = box.add_particle()
abox.x = (0,0,0)
abox.v = (0,0,0)
abox.q = (1,0,0,0)
abox.controlPoints.add_and_set_particles(xBF = points)
abox.RigidQuads.add_and_set_particles(vertexIndices = quads)

#------------------------------------------------------------------------------------------------------------------

mysim.run_until( endtime.get_value() )

#Verification that the output is identical to a reference value
from mpacts.tools.arrayhash import check_and_document_hash
check_and_document_hash( cells['x'], '6b1258f93e3cb6af2cac1b100963c73d', tolerance=1e-12)
#-----------------------------------------------------------------------------------------------------