Newer
Older
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
# Definitions applicable to many Spherocylinder models
# Individual models may need to override some of these definitions
# but sufficiently many models share these definitions that they are collected here
# spherocylinder models already have an explicit unitvector, even in 2D
unitvector(c::Spherocylinder) = c.backbone.eφ
mobilitytype(::Type{<:Spherocylinder}) = SpherocylinderMobility
"""
spawnchildren(parent::Spherocylinder, sim::Simulation)
Adds two new cells to the simulation, evenly spaced along the backbone of the parent
cell.
The new node separation is chosen to make the cells fit in length taken up by their parent
without.
The growth rate is randomly modified using the `growthrateretention` and a random
value drawn from the `growthratedistribution` using the simulation random number generator.
At the moment, all other parameters (except for the `id`) are inherited from the parent.
This method currently supports spherocylinders with and without the `attached` property.
If further properties need to be inherited, this method should be overridden.
(Or discuss, if it makes sense to modify this function definition)
"""
function spawnchildren(parent::T, sim::Simulation) where T <: Spherocylinder
# offset (in backbone parameter units) of the child cells
coffset = radius(parent)/(2*nodeseparation(parent))
# Note that forces are not continuous at division, but this doesn't seem to be that important
childlength = nodeseparation(parent)/2 - radius(parent)
for pos in (0.25 - coffset, 0.75 + coffset)
newcenter = parent.backbone(pos)
# do not spawn outside domain boundaries
isoutside(newcenter, sim.domain) && continue
# calculate new growth rate
γnew = parent.params.growthrateretention*parent.growthrate +
(1 - parent.params.growthrateretention)*rand(sim.rng, parent.params.growthratedistribution)
# optional kws in some models:
kws = Dict{Symbol, Any}()
hasproperty(parent, :attached) && (kws[:attached] = parent.attached)
# add child cell with new cell ID
addparticle!(sim, T;
centerpos = newcenter,
nodeseparation = childlength,
φ = parent.φ,
growthprog = 0.0,
growthrate = γnew,
_parent_id = parent.id,
params = parent.params,
kws...
)
end
end
"""
interactionrange(ptype::Type{Spherocylinder, sim::Simulation) → ir::Float64
Compute the range of the longest possible spatial interactions for [`Spherocylinder`](@ref)s.
Uses the largest `maxlength` in the registered parameter structs.
Spherocylinder models that feature e.g. an adhesion range should override this method.
"""
InPartS.interactionrange(::Type{<:Spherocylinder}, sim::Simulation) = maximum(p -> p.maxlength + 2p.radius, sim.params)
"""
spherocylinderforces!(cell1::Spherocylinder, cell2::Spherocylinder, sim::Simulation)
Computes the interaction forces between two cells using [`closestpoints()`](@ref).
Forces are distributed over `fP` and `fN` using [`distributeforce!()`](@ref).
"""
@forcedef function spherocylinderforces!(c1::Spherocylinder, c2::Spherocylinder, sim::Simulation)
# do simple square root free check first: Could they possibly touch if perfectly aligned?
δv = (c2.centerpos-c1.centerpos) % sim.domain
if δv⋅δv > (nodeseparation(c1)/2+radius(c1)+nodeseparation(c2)/2+radius(c2))^2
return
end
# find closest points modulo boundaries
s, t = closestpoints(c1.backbone, c2.backbone, sim.domain,
asquare = nodeseparation(c1) * nodeseparation(c1),
bsquare = nodeseparation(c2) * nodeseparation(c2))
# difference vector between closest points (again using MIC)
deltavector = (c1.backbone(s) - c2.backbone(t)) % sim.domain
# Avoid taking the square-root if particles are not actually touching
sqdistance = deltavector⋅deltavector
if sqdistance > (radius(c1)+radius(c2))^2 || sqdistance == 0
return
end
distance = sqrt(sqdistance)
# edge-to-edge distance
dedge = distance - radius(c1) - radius(c2)
# scalar force along unit vector
if dedge ≤ 0
# hertzian repulsion
force = hertzfactor(c1, c2, sim)*dedge*sqrt(-dedge)*deltavector/distance
else
# nothing more to do
return
end
# distribute forces over pseudo-nodes
distributeforce!(c1, -force, s)
distributeforce!(c2, force, t)
@force begin
midpoint = fold_back(c2.backbone(t) + deltavector ./2, sim.domain)
forcecontribution(c1.centerpos, -force, (midpoint-c1.centerpos)%sim.domain)
end
@force begin
midpoint = fold_back(c2.backbone(t) + deltavector ./2, sim.domain)
forcecontribution(c2.centerpos, force, (midpoint-c2.centerpos)%sim.domain)
end
return
end
"""
distributeforce!(c::Spherocylinder, force::SV{Float64}, r::Float64)
Distributes the given force over `fP` and `fN`, with `r*force` added to `fP` and
`(1-r)*force` added to `fN`. This is (of course) only applicable to spherocylinders
that implement node forces.
"""
@inline function distributeforce!(c::Spherocylinder{D}, force::SVector{D, Float64}, r) where {D}
c.fP += force * r
c.fN += force * (1 - r)
return
end
## Pretty printing
#one-line text output
Base.show(io::IO, cell::Spherocylinder) =
print(io, "$(typeof(cell))(r ≈ [$(join(nicestrlong.(cell.centerpos), ","))], l ≈ $(nicestrshort(nodeseparation(cell))), φ ≈ $(orientation_descriptor(cell)))")
# fancy multi-line output
Base.show(io::IO, ::MIME"text/plain", cell::Spherocylinder) = print(io,
"""$(typeof(cell)) $(repr(cell.id))
r ≈ [$(join(nicestrlong.(cell.centerpos), ","))]
l ≈ $(nicestrshort(nodeseparation(cell))); φ ≈ $(orientation_descriptor(cell))
g ≈ $(nicestrshort(cell.growthprog)); γ = $(nicestrshort(cell.growthrate))
$(hasproperty(cell, :attached) ? (cell.attached ? "attached\n" : "not attached\n") : "") \
params: $(cell.params)"""
)
function orientation_descriptor(c::AbstractRod)
if hasproperty(c, :φ)
return nicestrshort(c.φ)
elseif hasproperty(c, :backbone) && hasproperty(c.backbone, :eφ)
return "[$(join(nicestrlong.(c.backbone.eφ), ","))]"
else
return "unknown"
end
end