Source code for ecoli.library.lattice_utils

"""
==================================
Utilities for Lattice Environments
==================================
"""

import numpy as np
from scipy import constants
from vivarium.core.process import Process
from vivarium.library.units import remove_units
from vivarium.library.topology import get_in, assoc_path

from vivarium.library.units import units

AVOGADRO = constants.N_A / units.mol
# Cache units to save time
UNITS_UM = units.um
UNITS_MM = units.mM


[docs] def get_bin_site(location, n_bins, bounds): """Get a bin's indices in the lattice Parameters: location (list): A list of 2 floats that specify the x and y coordinates of a point inside the desired bin. n_bins (list): A list of 2 ints that specify the number of bins along the x and y axes, respectively. bounds (list): A list of 2 floats that define the dimensions of the lattice environment along the x and y axes, respectively. Returns: tuple: A 2-tuple of the x and y indices of the bin in the lattice. """ location = [coordinate.to(UNITS_UM).magnitude for coordinate in location] bounds = [bound.to(UNITS_UM).magnitude for bound in bounds] bin_site_no_rounding = np.array( [location[0] * n_bins[0] / bounds[0], location[1] * n_bins[1] / bounds[1]] ) bin_site = tuple(np.floor(bin_site_no_rounding).astype(int) % n_bins) return bin_site
[docs] def get_bin_volume(n_bins, bounds, depth): """Get a bin's volume Parameters: n_bins (list): A list of 2 ints that specify the number of bins along the x and y axes, respectively. bounds (list): A list of 2 floats that specify the lengths of the environment's sides in the x and y directions, respectively. In units of microns. depth (float): The depth of the environment, in microns. Returns: float: The volume of each bin in the lattice, in Liters. """ total_volume = depth * bounds[0] * bounds[1] return (total_volume / (n_bins[0] * n_bins[1])).to(units.L)
[docs] def count_to_concentration(count, bin_volume): """Convert a molecule count into a concentration. Parameters should all have units. Returned value will have units. Parameters: count (int): The number of molecules in the bin. bin_volume (float): The volume of the bin. Returns: float: The concentration of molecule in the bin. """ return count / (bin_volume * AVOGADRO)
[docs] def make_gradient(gradient, n_bins, size): """Create a gradient from a configuration **Random** A random gradient fills the field randomly with each molecule, with values between 0 and the concentrations specified. Example configuration: .. code-block:: python 'gradient': { 'type': 'random', 'molecules': { 'mol_id1': 1.0, 'mol_id2': 2.0 }}, **Uniform** A uniform gradient fills the field evenly with each molecule, at the concentrations specified. Example configuration: .. code-block:: python 'gradient': { 'type': 'uniform', 'molecules': { 'mol_id1': 1.0, 'mol_id2': 2.0 }}, **Gaussian** A gaussian gradient multiplies the base concentration of the given molecule by a gaussian function of distance from center and deviation. Distance is scaled by 1/1000 from microns to millimeters. Example configuration: .. code-block:: python 'gradient': { 'type': 'gaussian', 'molecules': { 'mol_id1':{ 'center': [0.25, 0.5], 'deviation': 30}, 'mol_id2': { 'center': [0.75, 0.5], 'deviation': 30} }}, **Linear** A linear gradient sets a site's concentration (c) of the given molecule as a function of distance (d) from center and slope (b), and base concentration (a). Distance is scaled by 1/1000 from microns to millimeters. .. math:: c = a + b * d Example configuration: .. code-block:: python 'gradient': { 'type': 'linear', 'molecules': { 'mol_id1':{ 'center': [0.0, 0.0], 'base': 0.1, 'slope': -10}, 'mol_id2': { 'center': [1.0, 1.0], 'base': 0.1, 'slope': -5} }}, **Exponential** An exponential gradient sets a site's concentration (c) of the given molecule as a function of distance (d) from center, with parameters base (b) and scale (a). Distance is scaled by 1/1000 from microns to millimeters. Note: base > 1 makes concentrations increase from the center. .. math:: c=a*b^d. Example configuration: .. code-block:: python 'gradient': { 'type': 'exponential', 'molecules': { 'mol_id1':{ 'center': [0.0, 0.0], 'base': 1+2e-4, 'scale': 1.0}, 'mol_id2': { 'center': [1.0, 1.0], 'base': 1+2e-4, 'scale' : 0.1} }}, Parameters: gradient: Configuration dictionary that includes the ``type`` key to specify the type of gradient to make. n_bins: A list of two elements that specify the number of bins to have along each axis. size: A list of two elements that specifies the size of the environment. """ bins_x = n_bins[0] bins_y = n_bins[1] length_x = size[0] length_y = size[1] fields = {} if gradient.get("type") == "random": for molecule_id, fill_value in gradient["molecules"].items(): field = fill_value * np.random.rand(bins_x, bins_y) fields[molecule_id] = field if gradient.get("type") == "gaussian": for molecule_id, specs in gradient["molecules"].items(): field = np.ones((bins_x, bins_y), dtype=np.float64) center = [specs["center"][0] * length_x, specs["center"][1] * length_y] deviation = specs["deviation"] for x_bin in range(bins_x): for y_bin in range(bins_y): # distance from middle of bin to center coordinates dx = (x_bin + 0.5) * length_x / bins_x - center[0] dy = (y_bin + 0.5) * length_y / bins_y - center[1] distance = np.sqrt(dx**2 + dy**2) scale = gaussian(deviation, (distance / 1000)) # multiply gradient by scale field[x_bin][y_bin] *= scale fields[molecule_id] = field elif gradient.get("type") == "linear": for molecule_id, specs in gradient["molecules"].items(): field = np.zeros((bins_x, bins_y), dtype=np.float64) center = [specs["center"][0] * length_x, specs["center"][1] * length_y] base = specs.get("base", 0.0) slope = specs["slope"] for x_bin in range(bins_x): for y_bin in range(bins_y): dx = (x_bin + 0.5) * length_x / bins_x - center[0] dy = (y_bin + 0.5) * length_y / bins_y - center[1] distance = np.sqrt(dx**2 + dy**2) field[x_bin][y_bin] += base + slope * (distance / 1000) fields[molecule_id] = field elif gradient.get("type") == "exponential": for molecule_id, specs in gradient["molecules"].items(): field = np.zeros((bins_x, bins_y), dtype=np.float64) center = [specs["center"][0] * length_x, specs["center"][1] * length_y] base = specs["base"] scale = specs.get("scale", 1) for x_bin in range(bins_x): for y_bin in range(bins_y): dx = (x_bin + 0.5) * length_x / bins_x - center[0] dy = (y_bin + 0.5) * length_y / bins_y - center[1] distance = np.sqrt(dx**2 + dy**2) field[x_bin][y_bin] = scale * base ** (distance / 1000) fields[molecule_id] = field elif gradient.get("type") == "uniform": for molecule_id, fill_value in gradient["molecules"].items(): fields[molecule_id] = np.full( (bins_x, bins_y), remove_units(fill_value), dtype=np.float64 ) return fields
[docs] def gaussian(deviation, distance): return np.exp(-np.power(distance, 2.0) / (2 * np.power(deviation, 2.0)))
[docs] def apply_exchanges( agents, fields, exchanges_path, location_path, n_bins, bounds, bin_volume ): # apply exchanges from each agent agent_updates = {} for agent_id, agent_state in agents.items(): exchanges = get_in(agent_state, exchanges_path) assert exchanges is not None location = get_in(agent_state, location_path) assert location is not None bin_site = get_bin_site(location, n_bins, bounds) reset_exchanges = {} for mol_id, value in exchanges.items(): # delta concentration exchange = value concentration = ( count_to_concentration(exchange, bin_volume).to(UNITS_MM).magnitude ) delta_field = np.zeros((n_bins[0], n_bins[1]), dtype=np.float64) delta_field[bin_site[0], bin_site[1]] += concentration fields[mol_id] += delta_field # reset the exchange value reset_exchanges[mol_id] = {"_value": -value, "_updater": "accumulate"} assoc_path(agent_updates, (agent_id,) + exchanges_path, reset_exchanges) return fields, agent_updates
[docs] class ExchangeAgent(Process): defaults = { "mol_ids": [], "default_exchange": 100, "max_move": 4, } def __init__(self, parameters=None): super().__init__(parameters)
[docs] def ports_schema(self): mol_ids = self.parameters["mol_ids"] def accumulate_location(value, update): return [ max(value[0] + update[0], 0 * units.um), max(value[1] + update[1], 0 * units.um), ] return { "boundary": { "external": {mol_id: {"_default": 0.0} for mol_id in mol_ids}, "exchanges": { mol_id: {"_default": self.parameters["default_exchange"]} for mol_id in mol_ids }, "location": { "_default": [0.5, 0.5], "_updater": accumulate_location, }, } }
[docs] def next_update(self, timestep, states): max_move = self.parameters["max_move"] mol_ids = self.parameters["mol_ids"] return { "boundary": { "exchanges": { mol_id: self.parameters["default_exchange"] for mol_id in mol_ids }, "location": [ np.random.uniform(-max_move, max_move) * units.um, np.random.uniform(-max_move, max_move) * units.um, ], } }
[docs] def make_diffusion_schema( exchanges_path, external_path, location_path, bounds, n_bins, depth, molecule_ids, default_field, ): # Place schemas at configured paths. agent_schema = {} location_schema = { "_default": [0.5 * bound for bound in bounds], "_emit": True, } assoc_path(agent_schema, location_path, location_schema) external_schema = { molecule: {"_default": 0.0 * units.mM} for molecule in molecule_ids } assoc_path(agent_schema, external_path, external_schema) exchanges_schema = {molecule: {"_default": 0.0} for molecule in molecule_ids} assoc_path(agent_schema, exchanges_path, exchanges_schema) # make the full schema schema = { "agents": {"*": agent_schema}, "fields": { field: { "_default": default_field.copy(), "_updater": "nonnegative_accumulate", "_emit": True, } for field in molecule_ids }, "dimensions": { "bounds": { "_default": bounds, "_updater": "set", "_emit": True, }, "n_bins": { "_default": n_bins, "_updater": "set", "_emit": True, }, "depth": { "_default": depth, "_updater": "set", "_emit": True, }, }, } return schema