Reference

The object of this class is the main API for running the simulation. It holds the state of a simulation and offers methods for advancing it.

Most attributes on this class describe the current state of molecules and cells. Both cells and molecules are referred to by index. Whenever cells are listed or represented in one dimension they are ordered consistently. E.g. at any point in the simulation world.cell_genomes[i] and world.cell_labels[i] refer to the genome and label of cell with index i. As cell number grows and shrinks cell indexes are not constant throughout the simulation. Whenever an operation modifies the number of cells (like kill_cells() or divide_cells()) cell indexes are updated. Likewise, molecules are always ordered as in Chemistry. E.g. world.chemistry.molecules[i] refers to the ith Molecule object and world.molecule_map[i] refers to a 2D map describing molecule concentrations of the ith molecule. Molecule indexes stay constant throughout a simulation.

The simulation works with lists and PyTorch Tensor representations of cell proteomes and molecule concentrations. You can get a Cell object from get_cell() with easy-to-interpret attributes (see Cell docs). However, for performance reasons only the lists and tensor representations are updated during a simulation. get_cell() is intended for analyzing simulation results.

For customizing a simulation one can interact with the working lists and tensors of the simulation. They could be used for reading cell and molecule information, or for modifying it. The following attributes are always updated during the simulation:

• cell_genomes: A list of strings representing cell genomes ordered by current cell index.
• labels: List of strings representing cell labels ordered by current cell index. Labels are strings that can be used to track cell origins. When spawning new cells (spawn_cell()) a random label is assigned to each cell. If a cell divides, its descendants get the same label.
• cell_map: Boolean 2D tensor describing which pixels on the map are occupied by a cell. Dimension 0 represents the x, dimension 1 the y direction.
• molecule_map: Float 3D tensor describing concentrations (in mM by default) for each Molecule. Dimension 0 represents the Molecule ordered as in Chemistry. Dimension 1 represents the x, dimension 2 the y direction. So, world.molecule_map[0, 1, 2] is the concentration of the 0th Molecule on pixel $(1, 2)$ of the world map.
• cell_molecules: Float 2D tensor describing intracellular concentrations (in mM by default) for each Molecule. Dimension 0 represents the cell ordered by current cell index. Dimension 1 represents the Molecule ordered as in Chemistry. So, world.cell_molecules[0, 1] represents the intracellular concentration in mM of the 1st Molecule in the cell which currently has index 0.
• cell_lifetimes: Integer 1D tensor describing how many time steps each cell survived since it was spawned or since its last division ordered by current cell index.
• cell_divisions: Integer 1D tensor describing how many times each cell's ancestors divided ordered by current cell index.

Some methods are provided to advance the simulation and interact with it:

• spawn_cells() spawn new cells from genomes and place them randomly on the map
• add_cells() add new cells from Cell objects and place them randomly on the map
• divide_cells() divide existing cells
• update_cells() update existing cells with a new genome and proteome
• kill_cells() kill existing cells
• move_cells() move existing cells to a random position in their Moore's neighborhood
• diffuse_molecules() let molecules diffuse and permeate by one time step
• degrade_molecules() let molecules degrade by one time step
• increment_cell_lifetimes() increment world.cell_lifetimes by 1
• enzymatic_activity() let cell proteins work for one time step
• get_cell() get a cell representation as a Cell object

Furthermore, there are methods for saving and loading a simulation. For any new simulation use save() once to save the whole world object to a pickle file (including its Chemistry, Genetics, Kinetics). You can restore it with from_file() later on. Then, to save the world's state you can use save_state(). This is a quick, lightweight save, but it only saves things that change during the simulation. Use load_state() to re-load a certain state.

PyTorch is heavily used in this simulation. Some important attributes are PyTorch Tensors (see above). By setting the device parameter most calculations can be moved to a GPU. This heavily increases performance and is recommended. When working with tensors one should try to use tensor methods. That way computations and data stay on the GPU. Here is a tutorial for working with tensors. Try to avoid sending data back and forth between CPU and GPU.

World object carries a Genetics instance on world.genetics, a Kinetics instance on world.kinetics, and a Chemistry instance on world.chemistry (which is just a reference to the Chemistry object that was used when initializing World). Usually, you don't need to touch them. But, if you want to override them or look into some details, see the docstrings of their classes for more information.

By default in this simulation molecule numbers can be thought of as being in mM, time steps in s, energies J/mol. Eventually, they are just numbers and can be interpreted as anything. However, with the default parameters in Kinetics and the way Molecule objects are defined it makes sense to interpret them in mM, s, and J.

Represents a molecule species which is part of the world, can diffuse, degrade, and be converted into other molecules.

Each molecule species which is supposed to be unique should have a unique name. In fact, if you initialize a molecule with the same name multiple times, only one instance of this molecule will be created. It allows you to define overlapping chemistries without creating multiple molecule instances of the same molecule species. You can also get a molecule instance from its name alone (Molecule.from_name()). However, this also means that if 2 molecules have the same name all attributes must match.

atp = Molecule("ATP", 10)
atp2 = Molecule("ATP", 10)
assert atp is atp2

atp3 = Molecule.from_name("ATP")
assert atp3 is atp

Molecule("ATP", 20)  # error because energy is different

By default in this simulation molecule numbers can be thought of as being in mM, time steps in s, energies in J/mol. Eventually, they are just numbers and can be interpreted as anything. However, together with the default parameters in Kinetics it makes sense to interpret them in mM, s, and J.

Molecule half life should represent the half life if the molecule is not actively deconstructed by a protein. Molecules degrade by one step whenever you call degrade_molecules(). You can setup the simulation to always call degrade_molecules() whenever a time step is finished.

Molecular diffusion in the 2D molecule map happens whenever you call diffuse_molecules(). The molecule map is world.molecule_map (see World). It is a 3D tensor where dimension 0 represents all molecule species of the simulation. They are ordered in the same way the attribute molecules is ordered in the Chemistry you defined. Dimension 1 represents x-positions and dimension 2 y-positions. Diffusion is implemented as a 2D convolution over the x-y tensor for each molecule species. This convolution has a 9x9 kernel. So, it alters the Moore's neighborhood of each pixel. How much of the center pixel's molecules are allowed to diffuse to the surrounding 8 pixels is defined by diffusivity. diffusivity is the ratio a/b when a is the amount of molecules diffusing to each of the 8 surrounding pixels, and b is the amount of molecules that stays on the center pixel. Thus, diffusivity=1.0 means all molecules of the center pixel are spread equally across the 9 pixels.

Molecules permeating cell membranes also happens with diffuse_molecules(). Cell molecules are defined in world.cell_molecules (see World). It is a 2D tensor where dimension 0 represents all cells and dimension 1 represents all molecule species. Again, molecule species are ordered in the same way the attribute molecules is ordered in the Chemistry you defined. Dimension 0 always changes its length depedning on how cells replicate or die. The cell index (cell.idx of Cell) for any cell equals the index in world.cell_molecules. E.g. the amount of molecule species currently in cell with index 100 is defined as world.cell_molecules[100]. permeability defines how much molecules can permeate from world.molecule_map into world.cell_molecules. Each cell lives on a certain pixel with x- and y-position. And although there are already molecules on this pixel, the cell has its own molecules. You could imagine the cell as a bag of molecule hovering over that pixel. permeability allows molecules from that pixel in the molecule map to permeate into the cell that lives on that pixel (and vice versa). So e.g. if cell 100 lives on pixel 12, 450 Molecules would be allowed to move from world.molecule_map[:, 12, 450] to world.cell_molecules[100, :]. Again, this happens separately for every molecule species depending on the permeability value. Specifically, permeability is the ratio of molecules that can permeate into the cell and the molecules that stay outside. Thus, a value of 1.0 means within one time step this molecule species spreads evenly between molecule map and cell.

Represents the chemistry with molecules and reactions available in a simulation.

molecules should include at least all molecule species that are mentioned in reactions. But it is possible to define more molecule species. Cells can use molecule species in transporer and regulatory domains, even if they are not included in any reaction. For stoichiometric coefficients > 1, list the molecule species multiple times. E.g. for 2A + B <-> C use reaction=([A, A, B], [C]) when A,B,C are molecule A, B, C instances.

Duplicate reactions and molecules will be removed on initialization. As any reaction can take place in both directions, it is not necessary to define both directions. You can use __and__ to combine multiple chemistries:

both = chemistry1 & chemistry2  # union of molecules and reactions

The chemistry object is used by World to know what molecule species exist. Reactions and molecule species are used to set up the world and create mappings for domains. On the world object there are some tensors that refer to molecule species (e.g. world.molecule_map). The molecule ordering in such tensors is always the same as the ordering in chemistry.molecules. E.g. if chemistry.molecules[2] is pyruvate, world.molecule_map[2] refers to pyruvate concentrations on the world molecule map. For convenience mappings are provided on this object:

• chemistry.mol_2_idx to map a Molecule object to its index
• chemistry.molname_2_idx to map a Molecule name string to its index

Protocol for domains

Object describing a catalytic domain.

The simulation works with tensor representations of cell proteomes. This class is a helper which makes interpreting a cell's proteome easier. You should not instantiate this class but instead get it from cell.proteome (see Cell).

Domain start and end describe the slice of the CDS python string. I.e. the index starts with 0, start is included, end is excluded. E.g. start=3, end=18 starts with the 4th and ends with the 18th base pair on the CDS.

Object describing a transporter domain.

The simulation works with tensor representations of cell proteomes. This class is a helper which makes interpreting a cell's proteome easier. You should not instantiate this class but instead get it from cell.proteome (see Cell).

is_exporter is only relevant in combination with other domains on the same protein. It defines in which transport direction this domain is energetically coupled with others.

Domain start and end describe the slice of the CDS python string. I.e. the index starts with 0, start is included, end is excluded. E.g. start=3, end=18 starts with the 4th and ends with the 18th base pair on the CDS.

Object describing a regulatory domain.

The simulation works with tensor representations of cell proteomes. This class is a helper which makes interpreting a cell's proteome easier. You should not instantiate this class but instead get it from cell.proteome (see Cell).

Domain start and end describe the slice of the CDS python string. I.e. the index starts with 0, start is included, end is excluded. E.g. start=3, end=18 starts with the 4th and ends with the 18th base pair on the CDS.

Object describing a protein.

The simulation works with tensor representations of cell proteomes. This class is a helper which makes interpreting a cell's proteome easier. You should not instantiate this class but instead get it from cell.proteome (see Cell).

CDS start and end describe the slice of the genome python string. I.e. the index starts with 0, start is included, end is excluded. E.g. cds_start=2, cds_end=31 starts with the 3rd and ends with the 31st base pair on the genome.

is_fwd describes whether the CDS is found on the forward (hypothetical 5'-3') or the reverse-complement (hypothetical 3'-5') side of the genome. cds_start and cds_end always describe the parsing direction / the direction of the hypothetical transcriptase. So, if you want to visualize a is_fwd=False CDS on the genome in 5'-3' direction you have to do n - cds_start and n - cds_stop if n is the genome length.

Object describing a cell and its environment.

The simulation works with tensor representations of cell proteomes and molecule concentrations. This class is a helper which makes interpreting a cell easier. You should not instantiate this class but instead get it from get_cell().

All parameters are exposed as attributes on the cell object. Some are calculated lazily just when they are accessed.

cell = world.get_cell(by_idx=3)
assert cell.idx == 5  # the cell's current index
cell.proteome  # proteome is calculated now

When a cell divides its genome and proteome are copied. Both descendants will recieve half of all molecules each. Both their n_divisions attributes are incremented. The cell's label will be copied as well. This way you can track how cells spread.

Protocol for domain factories

Factory for generating nucleotide sequences with CatalyticDomains.

For stoichiometric coefficients > 1, list the molecule species multiple times. E.g. for 2A + B <-> C use reaction=([A, A, B], [C]) when A,B,C are molecule A, B, C instances.

km and vmax are target values. Due to the way how codons are sampled for specific floats, the actual values for km and vmax might differ. The closest available value to the given value will be used.

If any optional argument is left out (None) it will be sampled randomly.

Factory for generating nucleotide sequences with TransporterDomains.

is_exporter is only relevant in combination with other domains on the same protein. It defines in which transport direction this domain is energetically coupled with others.

km and vmax are target values. Due to the way how codons are sampled for specific floats, the actual values for km and vmax might differ. The closest available value to the given value will be used.

If any optional argument is left out (None) it will be sampled randomly.

Factory for generating nucleotide sequences with RegulatoryDomains.

km is a target value. Due to the way how codons are sampled for specific values, the final value for km might differ. The closest available value to the given value will be used.

If any optional argument is left out (None) it will be sampled randomly.

Factory for generating genomes that translate into a desired proteome.

The desired proteome is given as a list of proteins where each protein is represented by a list of domain factories. As domain factories CatalyticDomainFact, TransporterDomainFact, RegulatoryDomainFact can be used. Use generate() to sample and generate a genome. There are always multiple base pair sequences which can encode the same proteome. Thus each generated genome might be different.

While the generated genome will always have the desired domains and proteins it might also include proteins which were not defined. This factory assembles a genome in one reading frame and direction. There is always the possibility that in another reading frame or on the reverse-complement of the genome other proteins are encoded. The larger the final genome, the higher the chances of this happening.

Class holding logic about transcribing and translating nucleotide sequence pairs.

During the simulation translate_genomes() is used. Its return value can be used to update cell parameters using Kinetics.set_cell_params() or to translate its abstract return value into a human readable Proteins description using Kinetics.get_proteome().

Translation happens only within coding sequences (CDSs). A CDS starts at every start codon and ends with the first in-frame encountered stop codon. Un-stopped CDSs are discarded. Both the forward and reverse complement of the nucleotide sequence are considered. Each CDS represents one Protein. All domains found within a CDS will be added as domains to that protein.

If you want to use your own genetics object for the simulation you can just assign it after creating World. E.g.:

world = World(chemistry=chemistry)
my_genetics = Genetics(p_transp_dom=0.1, stop_codons=("TGA",))
world.genetics = my_genetics

Class holding logic for simulating protein work. Usually this class is instantiated automatically when initializing World. You can access it on world.kinetics.

There are c cells, p proteins, s signals. Signals are basically Molecules, but we have to differentiate between intra- and extracellular molecules. So, there are twice as many signals as Molecules. Molecules are always ordered as in Chemistry: first, the intracellular ones, then the extracellular ones. Cells are ordered as in World. Proteins are ordered how they were found by Genetics.translate_genomes() in each cell.

Attributes on an object of this class describe cell parameters:

• Kmf, Kmb Affinities to all signals processed by each protein in each cell (c, p). Ratios of (f)orward and (b)ackward reaction affinities are their equilibrium constants.
• Kmr Affinity of each regulating signal for each protein and each cell (c, p, s) exponentiated by their hill coefficients.
• Vmax Maximum velocities of each protein in each cell (c, p).
• Ke Equilibrium constants derived from standard reaction free energy for each protein in each cell (c, p).
• N Stoichiometric number for each signal that is processed by each protein in each cell (c, p, s). Additionally, there are Nf and Nb which describe only forward and backward stoichiometric coefficients. This is needed in addition to N to properly describe reactions that involve co-factors (i.e. n=0).
• A Allosteric modulation for each signal in each protein in each cell (c, p, s). The number represents the hill coefficient. Positive coefficients have an activating effect, negative have an inhibiting effect.

The main method is integrate_signals(). When calling World.enzymatic_activity(), a matrix X of signals (c, s) is prepared and then integrate_signals(X) is called. Updated signals are returned and World.enzymatic_activity() writes them back to world.cell_molecules and world.molecule_map.

Another method, which ended up here, is set_cell_params() which reads proteomes and updates cell parameters accordingly. This is called whenever the genomes of some cells have changed. Currently, this is also the main performance bottleneck.

When this class is initialized it generates the mappings from indices to domain parameters by random sampling. These mappings are then used throughout the simulation. If you initialize this class again, these mappings will be different. Initializing World will also create one Kinetics instance on world.kinetics. If you want to access indices to domain mappings in your simulation, you should use world.kinetics.

Note: All default values are based on the assumption that energies are in J, a time step represents 1s, and molecule numbers are in mM. If you change the defaults, you need to reconsider how these numbers should be interpreted.

The kinetics used here can theoretically never create negative molecule concentrations or make a reaction overshoot its rquilibrium state. However, this simulation computes these things one step at a time. This and the numerical limits of data types used here can create edge cases that need to be dealth with: reaction quotients overshooting the equilibrium state and negative concentrations. Both are caused by high $V_{max}$ with low $K_m$ values. With the default $V_{max}$ and $K_m$ ranges and stoichiometric coefficients below 100, the heuristics in integrate_signals() are close enough for a meaningful simulation.

If an enzyme is far away from its equilibrium state $K_e$ and substrate concentrations are far above $K_m$ it will progress its reaction at full speed $V_{max}$. This rate can be so high that, within one step, $Q$ surpasses $K_e$. In the next step it will move at full speed into the opposite direction, overshooting $K_e$ again, and so on. Reactions with high stoichiometric coefficients are more prone to this as their rate functions are sharper. To combat this integrate_signals() works by iteratively approaching $V_{max}$ on multiple levels of the computation. The approach was tuned to compute fast and give satifisfying reults with certain conditions in mind. $V_{max} \le 100$, $K_m \ge 0.01$, and concentrations $X$ generally not much higher than 100. Violating these assumptions could lead to reaction quotients constantly overshooting their equilibrium.