Agar pad simulations

This notebook is currently being updated, and is unfinished!

[1]:
%load_ext autoreload
[2]:
%autoreload 2
[1]:
import numpy as np
import os
import pickle
from skimage.transform import rescale, rotate
import noise
import matplotlib.pyplot as plt
[2]:
import sys
sys.path.insert(1, '/home/georgeos/Documents/GitHub/SyMBac/') # Not needed if you installed SyMBac using pip
from SyMBac.drawing import raster_cell
from SyMBac.colony_simulation import ColonySimulation
/home/georgeos/Documents/GitHub/SyMBac/SyMBac/colony_simulation.py:8: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  from tqdm.autonotebook import tqdm

For microcolony simulations, we prefer to use CellModeller 1 instead of the inbuilt SyMBac simulator. We will be updating the following page with example models written for CellModeller, which cover simple microcolonies, and other microfluidic geometries: Example models.

1

Computational Modeling of Synthetic Microbial Biofilms Timothy J. Rudge, Paul J. Steiner, Andrew Phillips, and Jim Haseloff ACS Synthetic Biology 2012 1 (8), 345-352 DOI: 10.1021/sb300031n

[3]:
colonysim = ColonySimulation(
    cellmodeller_model= 'cellmodeller_ex1_simpleGrowth_modified.py',
    max_cells = 100,
    pix_mic_conv = 0.065,
    resize_amount = 3,
    save_dir = "test/",
)
[6]:
colonysim.run_cellmodeller_sim(num_sim=1)
Set up OpenCL context:
  Platform: NVIDIA CUDA
  Device: Quadro RTX 3000
Importing model cellmodeller_ex1_simpleGrowth_modified
      10           2 cells           0 contacts    0.000034 hour(s) or 0.002035 minute(s) or 0.122123 second(s)
      20           2 cells           0 contacts    0.000041 hour(s) or 0.002452 minute(s) or 0.147110 second(s)
      30           2 cells           0 contacts    0.000048 hour(s) or 0.002906 minute(s) or 0.174348 second(s)
      40           2 cells           0 contacts    0.000056 hour(s) or 0.003365 minute(s) or 0.201879 second(s)
      50           2 cells           0 contacts    0.000063 hour(s) or 0.003795 minute(s) or 0.227723 second(s)
      60           3 cells           1 contacts    0.000075 hour(s) or 0.004480 minute(s) or 0.268788 second(s)
   60     3 cells       2 cts       3 iterations  residual = 0.000513
      70           5 cells           4 contacts    0.000092 hour(s) or 0.005502 minute(s) or 0.330109 second(s)
   70     5 cells       4 cts       4 iterations  residual = 0.001992
      80           5 cells           4 contacts    0.000110 hour(s) or 0.006609 minute(s) or 0.396537 second(s)
   80     5 cells       4 cts       3 iterations  residual = 0.003237
      90           5 cells           5 contacts    0.000128 hour(s) or 0.007680 minute(s) or 0.460783 second(s)
   90     5 cells       4 cts       2 iterations  residual = 0.003866
     100           8 cells           7 contacts    0.000146 hour(s) or 0.008754 minute(s) or 0.525262 second(s)
  100     8 cells       7 cts       4 iterations  residual = 0.002115
     110           9 cells           8 contacts    0.000164 hour(s) or 0.009836 minute(s) or 0.590187 second(s)
  110     9 cells      10 cts       5 iterations  residual = 0.003253
     120          10 cells          14 contacts    0.000188 hour(s) or 0.011260 minute(s) or 0.675616 second(s)
  120    10 cells      14 cts       7 iterations  residual = 0.003114
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Input In [6], in <cell line: 1>()
----> 1 colonysim.run_cellmodeller_sim(num_sim=1)

File ~/Documents/GitHub/SyMBac/SyMBac/colony_simulation.py:56, in ColonySimulation.run_cellmodeller_sim(self, num_sim)
     54 # Run the simulation to ~n cells
     55 while len(self.simulation.cellStates) < self.max_cells:
---> 56     self.simulation.step()
     57 self.n_simulations += 1

File ~/miniconda3/envs/symbac/lib/python3.9/site-packages/CellModeller/Simulator.py:382, in Simulator.step(self)
    379             self.kill(state)
    381 self.phys.set_cells()
--> 382 while not self.phys.step(self.dt): #neighbours are current here
    383     pass
    384 if self.sig:

File ~/miniconda3/envs/symbac/lib/python3.9/site-packages/CellModeller/Biophysics/BacterialModels/CLBacterium.py:590, in CLBacterium.step(self, dt)
    588 if not self.progress_initialised:
    589     self.progress_init(dt)
--> 590 if self.progress():
    591     self.progress_finalise()
    592     return True

File ~/miniconda3/envs/symbac/lib/python3.9/site-packages/CellModeller/Biophysics/BacterialModels/CLBacterium.py:557, in CLBacterium.progress(self)
    555 def progress(self):
    556     if self.n_ticks:
--> 557         if self.tick(self.actual_dt):
    558             self.n_ticks -= 1
    559         return False

File ~/miniconda3/envs/symbac/lib/python3.9/site-packages/CellModeller/Biophysics/BacterialModels/CLBacterium.py:623, in CLBacterium.tick(self, dt)
    621 if not self.sub_tick_initialised:
    622     self.sub_tick_init(dt)
--> 623 if self.sub_tick(dt):
    624     self.sub_tick_finalise()
    625     return True

File ~/miniconda3/envs/symbac/lib/python3.9/site-packages/CellModeller/Biophysics/BacterialModels/CLBacterium.py:643, in CLBacterium.sub_tick(self, dt)
    641 self.build_matrix() # Calculate entries of the matrix
    642 #print "max cell contacts = %i"%cl_array.max(self.cell_n_cts_dev).get()
--> 643 self.CGSSolve(dt, alpha) # invert MTMx to find deltap
    644 self.add_impulse()
    645 return False

File ~/miniconda3/envs/symbac/lib/python3.9/site-packages/CellModeller/Biophysics/BacterialModels/CLBacterium.py:1030, in CLBacterium.CGSSolve(self, dt, alpha, substep)
   1026 max_iters = self.n_cells*7
   1028 for iter in range(max_iters):
   1029     # Ap
-> 1030     self.calculate_Ax(self.Ap_dev[0:self.n_cells], self.p_dev[0:self.n_cells], dt, alpha)
   1032     # p^TAp
   1033     pAp = self.vdot(self.p_dev[0:self.n_cells], self.Ap_dev[0:self.n_cells]).get()

File ~/miniconda3/envs/symbac/lib/python3.9/site-packages/CellModeller/Biophysics/BacterialModels/CLBacterium.py:951, in CLBacterium.calculate_Ax(self, Ax, x, dt, alpha)
    939 def calculate_Ax(self, Ax, x, dt, alpha):
    941     self.program.calculate_Bx(self.queue,
    942                               (self.n_cells, self.max_contacts),
    943                               None,
   (...)
    949                               x.data,
    950                               self.Mx_dev.data).wait()
--> 951     self.program.calculate_BTBx(self.queue,
    952                                 (self.n_cells,),
    953                                 None,
    954                                 numpy.int32(self.max_contacts),
    955                                 self.cell_n_cts_dev.data,
    956                                 self.n_cell_tos_dev.data,
    957                                 self.cell_tos_dev.data,
    958                                 self.fr_ents_dev.data,
    959                                 self.to_ents_dev.data,
    960                                 self.Mx_dev.data,
    961                                 Ax.data).wait()
    962     # Tikhonov test
    963     #self.vaddkx(Ax, numpy.float32(0.01), Ax, x)
    964
    965     # Energy mimizing regularization
    966     self.program.calculate_Mx(self.queue,
    967                                   (self.n_cells,),
    968                                   None,
   (...)
    974                                   x.data,
    975                                   self.Mx_dev.data).wait()

File ~/miniconda3/envs/symbac/lib/python3.9/site-packages/pyopencl/__init__.py:901, in _add_functionality.<locals>.kernel_call(self, queue, global_size, local_size, *args, **kwargs)
    895 def kernel_call(self, queue, global_size, local_size, *args, **kwargs):
    896     # __call__ can't be overridden directly, so we need this
    897     # trampoline hack.
    898
    899     # Note: This is only used for the generic __call__, before
    900     # kernel_set_scalar_arg_dtypes is called.
--> 901     return self._enqueue(self, queue, global_size, local_size, *args, **kwargs)

File <pyopencl invoker for 'calculate_BTBx'>:8, in enqueue_knl_calculate_BTBx(self, queue, global_size, local_size, arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, global_offset, g_times_l, allow_empty_ndrange, wait_for)

KeyboardInterrupt:
[ ]:
colonysim.get_simulation_dirs()
[ ]:
pickles = colonysim.get_simulation_pickles()
[ ]:
colonysim.get_max_scene_size()
[ ]:
colonysim.draw_simulation_OPL(n_jobs = -1, FL=True, density = 0.1, random_distribution = "uniform", distribution_args = (0.9, 3))
[ ]:
from SyMBac.colony_renderer import ColonyRenderer
[ ]:
from SyMBac.PSF import PSF_generator
from SyMBac.renderer import convolve_rescale
from skimage.util import random_noise
from scipy.ndimage import gaussian_filter

from skimage.exposure import rescale_intensity

[ ]:
my_kernel = PSF_generator(
    radius = 50,
    wavelength = 0.75,
    NA = 1.45,
    n = 1.4,
    resize_amount = 3,
    pix_mic_conv = 0.065,
    apo_sigma = 8,
    mode="simple fluo",
    condenser = "Ph3",
    offset = 0.02
)
my_kernel.calculate_PSF()
my_kernel.plot_PSF()

[ ]:
my_renderer = ColonyRenderer(colonysim, my_kernel)
[ ]:
test_img = my_renderer.render_scene(-1)
[ ]:
mask = my_renderer.mask_loader(-1)
[ ]:
plt.imshow(test_img, cmap="Greys_r")
[ ]:
plt.imshow(np.roll(test_img, [0,], axis=[1]), cmap="Greys_r")
[ ]:
my_renderer.generate_random_samples(1000, 0.2, "training_data")
[ ]:
import random
random.choice([(0, mask.shape[0]), (1, mask.shape[1]), ([0,1], mask.shape)])
[ ]: