Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit dc5e675c by Jigyasa Watwani
Options

more flags

parent edf05acf
Inline Side-by-side
Showing with 164 additions and 60 deletions
growing_domain/full_model/growing_domain.py
...@@ -56,9 +56,20 @@ class Growing_Domain(object): ...@@ -56,9 +56,20 @@ class Growing_Domain(object):
def advection(self, conc, vel, tconc): def advection(self, conc, vel, tconc):
return df.inner(df.div(vel * conc), tconc) return df.inner(df.div(vel * conc), tconc)
def active_stress(self, rho, c): def active_stress(self, rho, c, v):
if self.morphogen_mediated_active_stress: if self.morphogen_mediated_active_stress:
return self.lamda * rho/(rho + self.saturation_rho) * (c-self.average_c) * df.Identity(self.dimension) return self.lamda * rho/(rho + self.saturation_rho) * (c-self.average_c) * df.Identity(self.dimension)
if self.rho_t_dependent_active_stress:
if self.density_dependent_turnover_rho:
rho_t_dependence = (self.diffusion_rho*df.div(df.nabla_grad(rho))
- self.turnover_rho * rho * (rho - self.average_rho)
- df.div(rho * v))
else:
rho_t_dependence = (self.diffusion_rho*df.div(df.nabla_grad(rho))
- self.turnover_rho * (rho - self.average_rho)
- df.div(rho * v))
return self.lamda * rho_t_dependence/(rho + self.saturation_rho) * df.Identity(self.dimension)
else: else:
return self.lamda * rho/(rho + self.saturation_rho)* df.Identity(self.dimension) return self.lamda * rho/(rho + self.saturation_rho)* df.Identity(self.dimension)
...@@ -79,13 +90,15 @@ class Growing_Domain(object): ...@@ -79,13 +90,15 @@ class Growing_Domain(object):
return (elastic_stress + viscous_stress) return (elastic_stress + viscous_stress)
def stress(self, u, v, rho, c): def stress(self, u, v, rho, c):
return (self.passive_stress(u, v) + self.active_stress(rho, c)) return (self.passive_stress(u, v) + self.active_stress(rho, c, v))
def reaction_rho(self, rho, trho): def reaction_rho(self, rho, trho):
if self.morphogen_mediated_turnover_rho: if self.morphogen_mediated_turnover_rho:
return self.average_c * df.inner(rho - self.average_rho, trho) return self.average_c * df.inner(rho - self.average_rho, trho)
elif self.morphogen_mediated_average_rho: elif self.morphogen_mediated_average_rho:
return self.turnover_rho * df.inner(rho - self.average_c, trho) return self.turnover_rho * df.inner(rho - self.average_c, trho)
elif self.density_dependent_turnover_rho: # to model that cells are born or die only where the density is non-zero
return self.turnover_rho * df.inner(rho * (rho - self.average_rho), trho)
else: else:
return self.turnover_rho * df.inner(rho - self.average_rho, trho) return self.turnover_rho * df.inner(rho - self.average_rho, trho)
...@@ -115,21 +128,21 @@ class Growing_Domain(object): ...@@ -115,21 +128,21 @@ class Growing_Domain(object):
r2 = "+".join(['x['+str(i)+']*x['+str(i)+']' for i in range(self.dimension)]) r2 = "+".join(['x['+str(i)+']*x['+str(i)+']' for i in range(self.dimension)])
if self.dimension==1: if self.dimension==1:
rho0='1-cos(2*pi*x[0]/L)' rho0='1-cos(2*pi*x[0]/R)'
else: else:
# circular domain --> circular domain # circular domain --> circular domain
# rho0 = '1 + cos(pi*%s/R)' %r2 rho0 = '(a/2)*(1-tanh((sqrt(%s) - r0)/w))' %r2
rho0 = 'a*(1-tanh((sqrt(%s) - r0)/w))' %r2 # rho0 = 'a*(1-tanh((sqrt(%s) - r0)/w))' %r2 # grows to a different size
# circular domain --> elliptical domain (w=0.3, lambda = -30, k=1) # circular domain --> elliptical domain (w=0.19, lambda = -10, k=1)
# rho0 = 'a*x[0]*x[0]*(1-tanh((sqrt(%s) - r0)/w))' %r2 # rho0 = 'a*x[0]*x[0]*(1-tanh((sqrt(%s) - r0)/w))' %r2
# circular domain --> dumbbell shaped domain (w = 0.8, k = 0.01, lambda = -80) # circular domain --> dumbbell shaped domain (w = 0.8, k = 0.01, lambda = -80)
# rho0 = '0.25 + a*x[0]*x[1]*(1-tanh((sqrt(%s) - r0)/w))' %r2 # rho0 = '0.25 + a*x[0]*x[1]*(1-tanh((sqrt(%s) - r0)/w))' %r2
# circular --> kidney (w=0.8, r0=0.5, k=0.1, lambda=-30) # circular --> kidney (w=0.8, r0=0.5, k=0.1, lambda=-30)
# rho0 = '0.6 + x[0]*(1-tanh((sqrt(%s) - r0)/w))' %r2 # rho0 = 'x[0]*(1-tanh((sqrt(%s) - r0)/w))' %r2
rho0 = df.interpolate(df.Expression(rho0, pi=np.pi, R=self.system_size, a=self.average_rho, r0=0.5, w=0.3, degree=1), self.function_space.sub(2).collapse()) rho0 = df.interpolate(df.Expression(rho0, pi=np.pi, R=self.system_size, a=self.average_rho, r0=0.5, w=0.3, degree=1), self.function_space.sub(2).collapse())
# add noise # add noise
...@@ -138,9 +151,11 @@ class Growing_Domain(object): ...@@ -138,9 +151,11 @@ class Growing_Domain(object):
# ic for c # ic for c
if self.morphogen: if self.morphogen:
c0 = '1 + 0.1*a*cos(2.0*pi*%s/R)' %r2 # c0 = '1 + 0.1*a*cos(2.0*pi*sqrt(%s)/R)' %r2 # shrinks and stops
# c0 = df.interpolate(df.Constant(0.0), self.function_space.sub(3).collapse()) c0 = '1 + 0.1*a*cos(2.0*pi*%s/R)' %r2 # grows and stops
c0 = df.interpolate(df.Expression(c0,a = self.average_c, R=self.system_size, pi=np.pi, degree=1), self.function_space.sub(3).collapse()) # c0 = '1 + 0.1*a*cos(2.0*pi*x[0]/R)' # grows very slightly, shrinks even more slightly, changes shape
c0 = df.interpolate(df.Expression(c0, a = self.average_c, R=self.system_size, pi=np.pi, degree=1), self.function_space.sub(3).collapse())
# add noise # add noise
noise_c = self.noise_level * (2*np.random.random(c0.vector().size())-1) noise_c = self.noise_level * (2*np.random.random(c0.vector().size())-1)
c0.vector()[:] = c0.vector()[:] + noise_c c0.vector()[:] = c0.vector()[:] + noise_c
...@@ -219,9 +234,6 @@ class Growing_Domain(object): ...@@ -219,9 +234,6 @@ class Growing_Domain(object):
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True) self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True) self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
radius_array = np.zeros(maxsteps+1)
radius_array[0] = np.max(np.sqrt(self.mesh.coordinates()[:][0]**2 + self.mesh.coordinates()[:][1]**2))
for steps in progressbar.progressbar(range(1, maxsteps + 1)): for steps in progressbar.progressbar(range(1, maxsteps + 1)):
# solve # solve
df.solve(self.form == 0, self.f, self.bc) df.solve(self.form == 0, self.f, self.bc)
...@@ -243,7 +255,6 @@ class Growing_Domain(object): ...@@ -243,7 +255,6 @@ class Growing_Domain(object):
# move mesh # move mesh
dr = df.project(u, self.function_space.sub(0).collapse()) dr = df.project(u, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr) df.ALE.move(self.mesh, dr)
radius_array[steps] = np.max(np.sqrt(self.mesh.coordinates()[:][0]**2 + self.mesh.coordinates()[:][1]**2))
# update time # update time
self.time += self.timestep self.time += self.timestep
self.uFile.close() self.uFile.close()
...@@ -252,7 +263,6 @@ class Growing_Domain(object): ...@@ -252,7 +263,6 @@ class Growing_Domain(object):
if self.morphogen: if self.morphogen:
self.cFile.close() self.cFile.close()
return (radius_array)
if __name__ == '__main__': if __name__ == '__main__':
import json, datetime import json, datetime
...@@ -270,15 +280,7 @@ if __name__ == '__main__': ...@@ -270,15 +280,7 @@ if __name__ == '__main__':
params['timestamp'] = timestamp params['timestamp'] = timestamp
gd = Growing_Domain(params) gd = Growing_Domain(params)
gd.solve()
radius_array = gd.solve()
times = np.linspace(0, params['maxtime'], int(params['maxtime']/params['timestep']) + 1 )
fig,ax = plt.subplots(1,1)
ax.set_yscale("log")
ax.set_ylabel("$log R(t)$")
ax.set_xlabel("$t$")
ax.scatter(times, radius_array, facecolors='none', edgecolors='b')
plt.show()
with open(params['timestamp'] + '_parameters.json', "w") as fp: with open(params['timestamp'] + '_parameters.json', "w") as fp:
json.dump(params, fp, indent=4) json.dump(params, fp, indent=4)
......
growing_domain/full_model/parameters.json 0 → 100644
{
"morphogen" : false,
"dimension" : 2,
"resolution" : 100,
"system_size" : 1,
"meshfile" : "disk.xml.gz",
"morphogen_mediated_active_stress": false,
"rho_t_dependent_active_stress" : true,
"bulk_elasticity" : 1.0,
"shear_elasticity": 1.0,
"shear_viscosity" : 1.0,
"bulk_viscosity" : 1.0,
"friction" : 1.0,
"lamda": -0.1,
"diffusion_rho" : 1.0,
"morphogen_mediated_turnover_rho": false,
"density_dependent_turnover_rho": false,
"morphogen_mediated_average_rho": false,
"turnover_rho" : 1.0,
"average_rho" : 1.0,
"saturation_rho" : 1.0,
"diffusion_c" : 1.0,
"turnover_c" : 1.0,
"average_c" : 1.0,
"noise_level": 0.0,
"timestep" : 0.1,
"savetime": 0.2,
"maxtime" : 40
}
\ No newline at end of file
growing_domain/full_model/viz_growing_domain.py
...@@ -83,16 +83,19 @@ def get_data(params, DIR=''): ...@@ -83,16 +83,19 @@ def get_data(params, DIR=''):
def visualize(params, DIR='', offscreen=False): def visualize(params, DIR='', offscreen=False):
if params['morphogen']: if params['morphogen']:
(times, topology, geometry, u, v, rho, c) = get_data(params, DIR) (times, topology, geometry, u, v, rho, c) = get_data(params, DIR)
else: else:
(times, topology, geometry, u, v, rho) = get_data(params, DIR) (times, topology, geometry, u, v, rho) = get_data(params, DIR)
radius = np.linalg.norm(geometry, axis=2)
radial_coordinate = np.linalg.norm(geometry, axis=2)
if params['dimension']==1: if params['dimension']==1:
# c(x,t) and rho(x,t) vs x and t
if params['morphogen']: if params['morphogen']:
fig, (axrho, axc) = plt.subplots(2,1, sharex=True, figsize=(8,8)) fig, (axrho, axc) = plt.subplots(2,1, sharex=True, figsize=(8,8))
axc.set_xlabel(r'$x$') axc.set_xlabel(r'$x$')
axc.set_ylabel(r"$c(x,t)$") axc.set_ylabel(r"$c(x,t)$")
axc.set_xlim(np.min(radius), np.max(radius)) axc.set_xlim(np.min(radial_coordinate), np.max(radial_coordinate))
axc.set_ylim(np.min(c), np.max(c)) axc.set_ylim(np.min(c), np.max(c))
else: else:
...@@ -100,20 +103,20 @@ def visualize(params, DIR='', offscreen=False): ...@@ -100,20 +103,20 @@ def visualize(params, DIR='', offscreen=False):
axrho.set_xlabel(r'$x$') axrho.set_xlabel(r'$x$')
axrho.set_ylabel(r"$\rho(x,t)$") axrho.set_ylabel(r"$\rho(x,t)$")
axrho.set_xlim(np.min(radius), np.max(radius)) axrho.set_xlim(np.min(radial_coordinate), np.max(radial_coordinate))
axrho.set_ylim(np.min(rho), np.max(rho)) axrho.set_ylim(np.min(rho), np.max(rho))
rhoplot, = axrho.plot(radius[0], rho[0],'o',ms=3) rhoplot, = axrho.plot(radial_coordinate[0], rho[0],'o',ms=3)
if params['morphogen']: if params['morphogen']:
cplot, = axc.plot(radius[0], c[0]) cplot, = axc.plot(radial_coordinate[0], c[0])
def update(value): def update(value):
ti = np.abs(times-value).argmin() ti = np.abs(times-value).argmin()
rhoplot.set_ydata(rho[ti]) rhoplot.set_ydata(rho[ti])
rhoplot.set_xdata(radius[ti]) rhoplot.set_xdata(radial_coordinate[ti])
if params['morphogen']: if params['morphogen']:
cplot.set_ydata(c[ti]) cplot.set_ydata(c[ti])
cplot.set_xdata(radius[ti]) cplot.set_xdata(radial_coordinate[ti])
plt.draw() plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02]) sax = plt.axes([0.1, 0.92, 0.7, 0.02])
...@@ -124,6 +127,20 @@ def visualize(params, DIR='', offscreen=False): ...@@ -124,6 +127,20 @@ def visualize(params, DIR='', offscreen=False):
slider.on_changed(update) slider.on_changed(update)
plt.show() plt.show()
# heat map
# L(t) vs t
length = np.array(np.zeros(len(times)))
for j in range(len(times)):
length[j] = np.max(radial_coordinate[j])
figlength, axlength = plt.subplots(1, 1,)
axlength.set_xlabel('$t$')
axlength.set_ylabel('$L(t)$')
axlength.set_xlim(np.min(times), np.max(times))
axlength.set_ylim(np.min(length), np.max(length))
axlength.plot(times, length)
plt.show()
else: else:
# heat map # heat map
n_cmap_vals = 16 n_cmap_vals = 16
...@@ -159,7 +176,25 @@ def visualize(params, DIR='', offscreen=False): ...@@ -159,7 +176,25 @@ def visualize(params, DIR='', offscreen=False):
xmin=times[0], xmax=times.max(), xmin=times[0], xmax=times.max(),
value=times[0], title=r"$t/\tau$") value=times[0], title=r"$t/\tau$")
vd.show(interactive=(not offscreen), zoom=0.8) vd.show(interactive=(not offscreen), zoom=0.2)
if offscreen:
FPS = 5
movFile = '%s_cell_density.mov' % params['timestamp']
fps = float(FPS)
command = "ffmpeg -y -r"
options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
tmp_dir = TemporaryDirectory()
get_filename = lambda x: os.path.join(tmp_dir.name, x)
for tt in range(len(times)):
idx = (abs(times-times[tt])).argmin()
update(idx)
slider.GetRepresentation().SetValue(times[tt])
fr = get_filename("%03d.png" % tt)
vd.io.screenshot(fr)
os.system(command + " " + str(fps)
+ " -i " + tmp_dir.name + os.sep
+ "%03d.png " + options + " " + movFile)
tmp_dir.cleanup()
if params['morphogen']: if params['morphogen']:
plotter1 = vd.plotter.Plotter(axes=0) plotter1 = vd.plotter.Plotter(axes=0)
...@@ -190,37 +225,67 @@ def visualize(params, DIR='', offscreen=False): ...@@ -190,37 +225,67 @@ def visualize(params, DIR='', offscreen=False):
vd.show(interactive=(not offscreen), zoom=0.8) vd.show(interactive=(not offscreen), zoom=0.8)
if offscreen: if offscreen:
FPS = 10 FPS1 = 5
movFile = '%s_numerical.mov' % params['timestamp'] movFile1 = '%s_morphogen.mov' % params['timestamp']
fps = float(FPS) fps1 = float(FPS1)
command = "ffmpeg -y -r" command1 = "ffmpeg -y -r"
options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4" options1 = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
tmp_dir = TemporaryDirectory() tmp_dir1 = TemporaryDirectory()
get_filename = lambda x: os.path.join(tmp_dir.name, x) get_filename1 = lambda x: os.path.join(tmp_dir1.name, x)
for tt in range(len(times)): for tt in range(len(times)):
idx = (abs(times-times[tt])).argmin() idx1 = (abs(times-times[tt])).argmin()
update(idx) update1(idx1)
slider.GetRepresentation().SetValue(times[tt]) slider1.GetRepresentation().SetValue(times[tt])
fr = get_filename("%03d.png" % tt) fr1 = get_filename1("%03d.png" % tt)
vd.io.screenshot(fr) vd.io.screenshot(fr1)
os.system(command + " " + str(fps) os.system(command1+ " " + str(fps1)
+ " -i " + tmp_dir.name + os.sep + " -i " + tmp_dir1.name + os.sep
+ "%03d.png " + options + " " + movFile) + "%03d.png " + options1 + " " + movFile1)
tmp_dir.cleanup() tmp_dir1.cleanup()
# R(t) vs t
radius = np.array(np.zeros(len(times)))
for j in range(len(times)):
radius[j] = np.max(radial_coordinate[j])
figradius1, axradius1 = plt.subplots(1,1)
axradius1.set_xlabel(r'$t$')
axradius1.set_xlim(np.min(times), np.max(times))
axradius1.set_ylim(np.min(radius), np.max(radius)+1)
axradius1.set_ylabel(r'$\log R(t)$')
axradius1.semilogy(times, radius)
figradius2, axradius2 = plt.subplots(1,1)
axradius2.set_xlabel(r'$t$')
axradius2.set_xlim(np.min(times), np.max(times))
# axradius2.set_ylim(np.min(radius), np.max(radius))
axradius2.set_ylabel(r'$R(t)$')
axradius2.plot(times, radius)
plt.show()
# c(r,t) vs r # c(r,t) vs r
fig1, axrho1 = plt.subplots(1,1, figsize=(8,8)) if params['morphogen']:
axrho1.set_xlabel(r'$r$') figradial, (axcradial, axrhoradial) = plt.subplots(2,1, figsize=(8,8))
axrho1.set_xlim(np.min(radius), np.max(radius)) axcradial.set_xlabel(r'$r$')
axrho1.set_ylim(np.min(rho), np.max(rho)) axcradial.set_xlim(np.min(radial_coordinate), np.max(radial_coordinate))
axrho1.set_ylabel(r'$c(r,t)$') axcradial.set_ylim(np.min(c), np.max(c))
cplot, = axrho1.plot(radius[0], rho[0],'o', ms=3) axcradial.set_ylabel(r'$c(r,t)$')
cplot, = axcradial.plot(radial_coordinate[0], c[0],'o', ms=3)
else:
figradial, axrhoradial = plt.subplots(1,1, figsize=(8,8))
axrhoradial.set_xlabel(r'$r$')
axrhoradial.set_xlim(np.min(radial_coordinate), np.max(radial_coordinate))
axrhoradial.set_ylim(np.min(rho), np.max(rho))
axrhoradial.set_ylabel(r'$\rho(r,t)$')
rhoplot, = axrhoradial.plot(radial_coordinate[0], rho[0], 'o', ms=3)
def update(value): def update(value):
ti = np.abs(times-value).argmin() ti = np.abs(times-value).argmin()
cplot.set_ydata(rho[ti]) rhoplot.set_ydata(rho[ti])
cplot.set_xdata(radius[ti]) rhoplot.set_xdata(radial_coordinate[ti])
if params['morphogen']:
cplot.set_ydata(c[ti])
cplot.set_xdata(radial_coordinate[ti])
plt.draw() plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02]) sax = plt.axes([0.1, 0.92, 0.7, 0.02])
...@@ -229,8 +294,8 @@ def visualize(params, DIR='', offscreen=False): ...@@ -229,8 +294,8 @@ def visualize(params, DIR='', offscreen=False):
fc='#999999') fc='#999999')
slider1.drawon = False slider1.drawon = False
slider1.on_changed(update) slider1.on_changed(update)
# plt.show() plt.show()
if __name__ == '__main__': if __name__ == '__main__':
import argparse, json import argparse, json
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment