cell_spreadingΒΆ

../_images/cell_spreading.png
'''
Mpacts script that simulates a deformable cell, defined as a cortical shell that dynamically adheres
to a rigid rectangular substrate.
'''
from mpacts.commands.onarrays.setvalue import SetValueCommand
from mpacts.commands.onarrays.transfer import SumArrayToParentCommand
from mpacts.contact.detectors.multigrid import MultiGridContactDetector
from mpacts.contact.detectors.fixedlist import FixedListContactDetector
from mpacts.contact.matrix.conjugategradient import DefaultConjugateGradientSolver
from mpacts.contact.models.bending.commonedge import ElasticShellBending, LinearShellBending
from mpacts.contact.models.collision.md.md_rt_matrix import MaugisDugdaleGeoFprimIntMatrix
from mpacts.contact.models.springs.elastic.elastic_matrix import ElasticShellTensionMatrix
from mpacts.core.arrays import create_array
from mpacts.core.simulation import Simulation
from mpacts.geometrygenerators.polyhedron import Connectivity
from mpacts.geometrygenerators.trianglegeometries import UnitIcoSphere
from mpacts.core.valueproperties import Variable, VariableFunction
import mpacts.commands.force.constraints.areaconservation as area
import mpacts.commands.force.constraints.volumeconservation as vol
import mpacts.contact.matrix.pythontools as cmpt
import mpacts.geometrygenerators.quadgeometries as quadgeo
from mpacts.tools.setcontactdatastorage import SetContactDataStorage
from mpacts.core.units import unit_registry as u
import DEMutilities.postprocessing.h5storage as h5s
import mpacts.commands.monitors.progress as pp
from mpacts.io.datasave import DataSaveCommand
import mpacts.particles as prt
import numpy as np
#----------------------------------------------------------------------------------------------------------------

#We start by making the 'root', a simulation object, from which all simulation branches will emerge
mysim = Simulation("simulation", timestep=1)
SetContactDataStorage( "Vector" )#Lay-out of contact data in memory.
params = mysim('params')
#!/usr/bin/env python-mpacts

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

#Defining variables (geometrical/interaction/material/numerical properties):
#The 'natural' units of the cell system are 'kPa', 'um' (micrometer), 'nN' (nanonewton) and s (seconds)
#NOTE these are not realistic parameters for biological cells, but are quick enough to simulate and are in the
#proper order of magnitude

#Mechanics
e_cortex   = Variable( "e_cortex"        , params, value=5*u('kPa') )#Young's modulus of the cortex
ec_cortex  = Variable( "ec_cortex"       , params, value=50*u('kPa'))#Contact stiffness (really high)
nu_cortex  = Variable( "nu_cortex"       , params, value=1./3. )#Poisson number of the cortex
kv_cell    = Variable( "kv_cell"         , params, value=0.5*u('kPa') )#Effective bulk modulus of the cell
ka_cortex  = Variable( "ka_cortex"       , params, value=0.2*u('nN/um') )#Local area conservation modulus
kd_cortex  = Variable( "kd_cortex"       , params, value=0.2*u('nN/um') )#Global area conservation modulus
mu_cortex  = Variable( "mu_cortex"       , params, value=0.5*u('kPa*s   '))#Viscosity of the cortex
t_cortex   = Variable( "t_cortex"        , params, value=1.0*u('um'))#Thickness of the cortex
visc_liq   = Variable( "mu_liquid"       , params, value=0.04*u('kPa*s') )#Effective viscosity of the cytosol
st_cortex  = Variable( "st_cortex"       , params, value=0.25*u('nN/um'))#Active surface tension (contractility)

e_plane    = Variable( "e_plane"         , params, value=1000*u('kPa'))#Young's modulus of the substrate
nu_plane   = Variable( "nu_plane"        , params, value=0.5 )#Poisson number of the substrate

gn_contact = Variable( "gamma_n_contact", params, value=2.0*u('kPa*s/um') )#Normal contact friction
gt_contact = Variable( "gamma_t_contact", params, value=1.0*u('kPa*s/um') )#Tangential (sliding) contact friction
adhesion   = Variable( "adhesion"       , params, value=0.5*u('nN/um') )#Cell-substrate adhesion energy
erange     = Variable( "effective_range", params, value=0.1*u('um') )#Effective adhesive range in MD-model
d_init     = Variable( 'initial_overlap', params, value=0.03*u('um') )#Initial overlap of cell with substrate

#Geometry
radius     = Variable( "radius_cell"     , params, value=6*u('um') )#Radius of the cell
subd       = Variable( "subdivision"     , params, value=4 )#Degree of mesh refinement as the icosahedron subdivision
lplane     = Variable( "l_plane"         , params, value=25*u('um') )#Side length of the plane

#Simulation
dt         = Variable( "timestep"        , params, value=2.5*u('ms'))#Simulation time step
endtime    = Variable( "endtime"         , params, value=6.*u('s') )#Tota l time to simulate
cd_ue      = Variable( "cd_update_every" , params, value=25 )#Update verlet lists every 'cd_update_every' steps
cd_kd      = Variable( "cd_keep_distance", params, value=1.0*u('um') )#margin to keep contacts for contact detection
cg_tol     = Variable( "cg_tolerance"    , params, value=1e-4)#Conjugate gradient rhs residual (force)
out_int    = Variable( "output_interval" , params, value=150*u('ms'))#Write an output file at this rate
rhv        = Variable( "rhv"             , params, value=1e8  )#'Really High Value' where needed

gn_cortex  = VariableFunction( "gamma_n_cortex"  , params, function='$t_cortex$*$mu_cortex$')
gt_cortex  = VariableFunction( "gamma_t_cortex"  , params, function='$t_cortex$*$mu_cortex$')
gamma_liq  = VariableFunction( "gamma_liquid", params, function ='1.5*$mu_liquid$/$radius_cell$' )
rhvmat     = VariableFunction( "gamma_plate", params, function='( $rhv$,0,$rhv$,0,0,$rhv$ )')
max_r      = VariableFunction( "max_inv_curv", params, function='10*$radius_cell$' )

#We may not forget to actually set the simulation timestep with a variable
mysim.set(timestep = dt)

#-----------------------------------------------------------------------------------------------------------------
#Some calculations with our parameters:
storage = h5s.H5Storage(mysim.name()+'.h5', openmode='wa')
storage.add_simulation_info( mysim )

#------------------------------------------------------------------------------------------------------------------
#Setting up the particle container lay-out
cells = prt.ParticleContainer("cell"
                             , prt.DeformableBody.compose(  ( prt.Node,"nodes" )
                                                         ,  ( prt.DeformableRoundedTriangle
                                                            , "triangles"
                                                            , "nodes" )
                                                         )
                             , mysim )

ftriangle = create_array("Vector", "F", cells('triangles'))
contact_area = create_array("Scalar","contact_area",cells('triangles'))
ca_parent = create_array("Scalar","contact_area", cells)
x_parent = create_array("Point", "x", cells)
p0 = create_array("Scalar", "internal_pressure", cells )
create_array("Scalar", "thickness", cells("triangles") )

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

create_array("Scalar", "contact_area", plane("RigidQuads"))
create_array("Vector", "F", plane("RigidQuads"))

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

#Adding 'physics': commands that are executed everytime, that add mechanical behavior/geometrical calculations.

SetValueCommand("ZeroCellContactArea",mysim, value=0, array=ca_parent)
SetValueCommand("ZeroCellPosition", mysim, value=(0,0,0), array=x_parent)
SetValueCommand("InitializeShellThickness", mysim, value=t_cortex, array=cells('triangles')['thickness'])

cells("triangles").SumArrayToParentCmd( name="SumContactAreaToParentCmd"
                                      , array=contact_area, array_parent = ca_parent )
cells('nodes').SumArrayToParentCmd( name = "SumNodePositionsToParentCmd"
                                  , array = cells('nodes')['x']
                                  , array_parent = x_parent )
cells("nodes").CountPerParentCmd()
cells('nodes').NormalizeParentArrayWithChildCountCmd( name = "DivideCellPositionWithNodeCountCmd"
                                                    , array = cells('x')
                                                    , array_count = cells['count_nodes'] )

SetValueCommand("ZeroTriangleContactArea",mysim, value=0, array=contact_area)
SetValueCommand("ZeroTriangleForce",mysim,value=(0,0,0),array=ftriangle)
cells("nodes").SetContactMatrixDiagonalCmd( visc=0 )#We will add to this later, set it to zero first!

cells("triangles").TriangleNormalsAndAreaCmd(x=cells('nodes')['x'])
cells("triangles").ComputeCurvatureCmd( maximal_radius=max_r,enforce_positive=True, save_for_nodes=True )
area.AreaConservation(  cells, cells("triangles"), ka=ka_cortex, kd=kd_cortex )
vol.VolumeConservation( cells, cells("triangles"), kv=kv_cell )

cells("triangles").TriangleCentersCmd( x=cells('nodes')['x'] )
cells("triangles").EncompassingSpheresCmd(x=cells('nodes')['x'])

cells("nodes").TimeIntegration_ForwardEuler_Generic_Cmd()
cells("nodes").FrictionOnAreaCmd( area = cells('nodes')['area']
                                , gamma_normal = gamma_liq
                                , gamma_tangential = gamma_liq )

cells("nodes").ComputeSurfaceTensionCmd( surface_tension = st_cortex
                                       , base_pressure = cells['internal_pressure'])

SetValueCommand("ZeroTensionNodes", mysim, value=0, array=cells("nodes")['tension'])

cells("triangles").VTKWriterCmd( executeInterval=out_int, select_all = True)
cells.CSVWriterCmd( executeInterval=out_int )
#cells('triangles').CSVWriterCmd(executeInterval=out_int)//Use this to visualize e.g. encompassing spheres/triangle

plane("RigidQuads").QuadNormalsCmd( x=plane('controlPoints')['x'] )
plane("RigidQuads").VTKWriterCmd( executeOnce=True, select_all = True )

SetValueCommand("ZerPlaneContactArea", mysim, value=0., array=plane("RigidQuads")['contact_area'])
SetValueCommand("ZerPlaneFprim", mysim, value=(0,0,0), array=plane("RigidQuads")['F'])
plane.SetContactMatrixDiagonalCmd( visc=rhvmat )

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

#Adding more 'physics' but this time as contact models: Interaction modules between different particles in simulation

spring = ElasticShellTensionMatrix( E = e_cortex
                                  , nu = nu_cortex
                                  , viscosity = mu_cortex
                                  , triangles = cells('triangles')
                                  , pc1 = cells('nodes')
                                  , pc2 = cells('nodes') )

bending = ElasticShellBending( E = e_cortex
                             , nu = nu_cortex
                             , pc1 = cells('triangles')
                             , pc2 = cells('triangles') )

CD_linear_spring = FixedListContactDetector( "CD_Springs", mysim, cmodel = spring )
CD_bending_triangles = FixedListContactDetector( "CD_Bending", mysim, cmodel = bending )

CM_cell_plane = MaugisDugdaleGeoFprimIntMatrix( pc1 = cells('triangles')
                                              , pc2 = plane("RigidQuads")
                                              , E1 = ec_cortex
                                              , E2 = e_plane
                                              , attrConst = adhesion
                                              , nu1= nu_cortex
                                              , nu2 = nu_plane
                                              , effective_range = erange
                                              , gamma_normal = gn_contact
                                              , gamma_tangential = gt_contact
                                              , Fprim1 = cells('triangles')['F']
                                              , Fprim2 = plane("RigidQuads")['F'])

CD_cell_plane = MultiGridContactDetector( "cell-plane", mysim
                                        , cmodel = CM_cell_plane
                                        , update_every = cd_ue
                                        , keep_distance = cd_kd )

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

ConjugateGradient = DefaultConjugateGradientSolver( mysim, tolerance = cg_tol, reset_x = False )

#Add a command that prints some information to the screen at a certain time interval:
printerlist = [ pp.DefaultProgressPrinter( mysim )
              , pp.PropertyPrinter( ConjugateGradient('steps'), "CG ")
              ]
printer = pp.PrinterChain( printerlist )
pp.ProgressIndicator("PrintProgress", mysim, printer, print_interval=5)

DataSaveCommand("DataSaveCommand", mysim, executeInterval = out_int )

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

#Actually adding the initial configuration of the cell, with nodes, triangles and connectivity ('springs')
print(" - subdividing icosphere")
p = UnitIcoSphere(int(subd.get_value()))
p.cleanLists()
print(" - calculating connectivity")
con = Connectivity(p)

xs = [(0,radius.get_value() - d_init.get_value(),0)]
rs = [radius.get_value()]

print(" - adding "+str(len(xs))+" cells")
for xsi,rsi in zip(xs,rs):
    cell = cells.add_particle()
    cell.internal_pressure = 2*st_cortex.get_value() / rsi#Makes sure that initial shape is equilibrium shape
    np_vertices = np.array(p.vertices)
    np_vertices *= rsi
    np_vertices += np.array(xsi)
    cell.nodes.add_and_set_particles( x=list(map(tuple, np_vertices)) )
    cell.triangles.add_and_set_particles( vertexIndices=p.triangles )
    cell.add_connectivity( CD_linear_spring, con.edgeCorners )
    cell.add_connectivity( CD_bending_triangles, con.edgeTriangles )


ctrlplane, vilplane = quadgeo.create_plane( lplane.get_value()
                                          , lplane.get_value()
                                          , perpendicular_to = 'y')
theplane = plane.add_particle()
theplane.x = (0,0,0)
theplane.q = (1,0,0,0)
theplane.controlPoints.add_and_set_particles( xBF = ctrlplane )
theplane.RigidQuads.add_and_set_particles( vertexIndices = vilplane )

#------------------------------------------------------------------------------------------------------------------
#All done now, we can start a valid simulation
print(" - starting simulation...")

print( mysim.print_command_tree() )
mysim.run_until(endtime.get_value())
storage.simulation_completed()
#------------------------------------------------------------------------------------------------------------------
#Verification
from mpacts.tools.arrayhash import check_and_document_hash
check_and_document_hash( cells('nodes')['x'], 'd1b575bdf2bdb745ab1c95bda73a114f', tolerance=1e-12)
#------------------------------------------------------------------------------------------------------------------