system_util

Short description

Compartment

class pyrid.system.system_util.Compartment(*args, **kwargs)

Compartment class contains attributes and methods used to describe a triangulated mesh compartment.

Attributes

triangle_idsint64[:]

triangle indices of the compartment mesh.

triangle_ids_exclTranspint64[:]

triangle indices of the compartment mesh, escluding triangles marked as transparent.

namestr

Name of the compartment.

AABBfloat64[:,:]

Axis Aligned Bounding Box of the compartment

originfloat64[:]

coordinates of AABB origin (AABB[0])

box_lengthsfloat64[:]

length of AABB in each dimension

centroidfloat64[:]

centroid of mesh

border_2darray_like 1-D

Information about edge of intersection between mesh and simualtion box

data-type: np.dtype([(‘edge_id’, np.int64), (‘triangle_ids’, np.int64),(‘direction_normal’, np.float64, (3,)), (‘edge_length’, np.float64),(‘cum_edge_length’, np.float64),], align=True)

border_lengthfloat64

length of edge of intersection between mesh and simualtion box

border_3darray_like 1-D

Information about faces of intersection between mesh and simulation box

data-type: np.dtype([(‘triangle_ids’, np.int64), (‘cum_area’, np.float64),], align=True)

border_areafloat64

area of intersection between mesh and simualtion box

groupsarray_like 1-D

triangle ids of face groups

data-type: np.dtype([(‘triangle_ids’, np.int64),], align=True)

areafloat64

Total area of mesh surface

volumefloat64

Total mesh volume

box_overlapbool

True is any intersection between mesh and simulation box

Methods

calc_centroid_radius_AABB(self, System):	calculates some general properties of the mesh compartment
add_border_2d(triangle_ids_border, System):	Registers the edge of intersection (if exists) between compartment mesh and simulation box.

add_border_2d(triangle_ids_border, System)

Registers the edge of intersection (if exists) between compartment mesh and simulation box. The intersection line is needed for fixed concdentration simulations as new molecules enter the simualtion box along this intersection line.

Note

See naming conventions for face groups!

Parameters

triangle_ids_borderint64[:]

triangle ids of the mesh triangle along the border

Systemobj

Instance of System class

add_border_3d(triangle_ids_border, System)

Registers the area of intersection (if exists) between compartment mesh and simulation box. The intersection area is needed for fixed concdentration simulations as new molecules enter the simualtion box along the faces of the intersection.

Note

See naming conventions for face groups!

Parameters

triangle_ids_borderint64[:]

triangle ids of the mesh triangle along the border

Systemobj

Instance of System class

add_group(group_name, triangle_ids_group, System)

Registers any triangle/face group that has been defined in the mesh ob file that is not a 2d or 3d border of intersectioj with the simulation box

Note

See naming conventions for face groups!

Parameters

group_namestr

name of the group

triangle_ids_groupint64[:]

triangle ids of the mesh triangle group

Systemobj

Instance of System class

calc_centroid_radius_AABB(System)

Calculates some general properties of the mesh compartment

• area

• volume

• centroid

• AABB (Axis Aligned Bounding Box), origin, box_lengths

The centroid (geometric center) of the mesh is calculated by averaging over all triangles centroids, weighted by their area.

Parameters

Systemobj

Instance of System class

MoleculeType

class pyrid.system.system_util.MoleculeType(*args, **kwargs)

A brief summary of the classes purpose and behavior

Attributes

namestr

Name of molecule

type_idint64

ID of molecule type

collision_typeint64 {0,1}

collision type of molecule. (Default = 0)

posfloat64[N,3]

positions of each particle

radiifloat64[N]

radii of particles

pos_rbfloat64[N,3]

position of each particle with exculded volume radius > 0.0

radii_rbfloat64[N]

radii of particles with exculded volume radius > 0.0

typesarray_like 1-D

array of names of particle types dtype : np.dtype(‘U20’)

radiusfloat64

Total radius of molecule

radius_2dfloat64

Total radius of molecule in plane (needed when distributing molecules on surface)

Dtransfloat64

Translational diffusion coefficient

Drotfloat64

Rotational diffusion coefficient

mu_rbfloat64[3,3]

Rotational mobility tensor

mu_tbfloat64[3,3]

Translational mobility tensor

mu_rb_sqrtfloat64[3,3]

Rotational mobility tensor

mu_tb_sqrtfloat64[3,3]

Square root of translational mobility tensor

mu_rb_sqrtfloat64[3,3]

Square root of rotational mobility tensor

D_ttfloat64[3,3]

Translational diffusion tensor

D_rrfloat64[3,3]

Rotational diffusion tensor

r_meanfloat64

Average radial displacement of diffusing molecule

volumefloat64

Volume of molecule

transition_rate_totalfloat64

Total transition rate of unimolecular reactions

um_reactionbool

True if unimolecular reaction has been defined for molecule type

number_umrint64

Number of unimolecular reactions

concentrationfloat64[:]

concentration of molecule in volume of each compartment outside simulation box (Needed for simulation with fixed concentration boundary condition)

concentration_surfacefloat64[:]

concentration of molecule on surface of each compartment outside simulation box (Needed for simulation with fixed concentration boundary)

l_perpfloat64

Mean diffusive displacement of molecule perpendicular to a plane.

h_membranefloat64

Distance of the molecules´ center of diffusion from the membrane surface.

Methods

update_um_reaction(rate)

Updates the the total transitaion rate for unimolecular reactions of this molecule type.

update_um_reaction(rate)

Updates the the total transitaion rate for unimolecular reactions of this molecule type.

Parameters

ratefloat64

Reaction rate

Notes

We cannot tell when a bimolecular reaction will occur, beacause we do not know when two diffusing molecules meet. However, for unimolecular reactions the expected lifetime distribution is known in case there is no interruption by any bimolecular reaction:

(29)\[\rho(t) = k*e^{-k t},\]

where \(k = \sum_i^n(k_i)\). As such, we can schedule a reaction event by using the upper equation Stiles and Bartol [11], Erban et al. [12]. For a unimolecular reaction occuring between t+dt we add this reaction to the reaction list evaluated at time t+dt. The reaction is then processed as any other bimolecular reaction that may occur within that timestep. If the unimolecular reaction is picked, it needs to be checked if different unimolecular reaction pathways exist. The fractional probablity for each of the n transitions is:

(30)\[p_i = k_i/\sum_n(k_j)\]

Based on the latter equation one of the n unimolecular transitions will be picked. As for bimolecular reactions, all other reactions that exist for the educt will be deleted. In case a bimolecular reaction is picked before the unimolecular reaction and in case it is successfull, the unimolecular reaction will be deleted from the reaction list and a new reaction time is drawn for the product molecule of the bimolecular reaction if any unimolecular reaction has been defined for the product type.

System

class pyrid.system.system_util.System(*args, **kwargs)

A brief summary of the classes purpose and behavior

Attributes

Nnb.int64

Total number of molecules in simulation

Npnb.int64

Total number of particles in simulation

Nmolnb.int64[:]

Number of molecules per molecule type

dimnb.int64 (3)

System dimension

dtnb.float64

Integration time step

box_lengthsnb.float64[:]

Length of simulation box in each dimension

volumenb.float64

Simulation box volume

Epotnb.float64[:]

Total potential energy

AABBnb.float64[2,3]

Axis Aligned Bounding Box of simulation box

AABB_allnb.float64[2,3]

Axis Aligned Bounding Box of simulation box plus mesh compartments (only different from simulation box if compartments overlap with box)

originnb.float64[3]

origin of simulation box (AABB[0])

emptynb.int64 (-1)

empty flag

verticesnb.float64[N,3]

Position of all vertices in each dimension

Mesharray_like

Mesh data

dtype : np.dtype([(‘triangles’, np.int64, (3,)), (‘edges’, np.int64, (3,2)), (‘border_edge’, np.int64, (3,)), (‘border_dim’, np.int64, (3,)), (‘border_normal’, np.float64, (3,3)), (‘neighbours’, np.int64, (3,)), (‘triangle_coord’, np.float64, (4,3)), (‘normal’, np.float64, (3,)), (‘triangle_distance’, np.float64), (‘triangle_centroid’, np.float64, (3,)), (‘triangle_area’, np.float64), (‘stamp’, np.uint64), (‘comp_id’, np.uint64), (‘group_id’, np.uint64),], align=True)

N_trianglesnb.int64

Total number of mesh triangles

triangle_idsnb.int64[:]

triangle ids of all compartments, excluding transparent faces (borders)

box_borderborder_box_lt

Triangle faces of the simulation box border (needed for distributing molecules in case of fixed concentration boundary condition)

box_areanb.float64

Area of triangle faces of the simulation box border

boundary_2dnb.types.boolean

True if intersection of compartment with simulation box exists.

triangle_ids_transparentnb.int64[:]

Triangle ids of transparent faces (particles do not collide)

Compartmentsnb.types.DictType(Comp_keytype, Comp_valuetype)

Dictionary containing the Compartment instances

n_compsnb.int64

Total number of compartments

areanb.float64

Total surface area of all compartments

cells_per_dimnb.int64[3]

Number of cells in each dimension of system grid

cell_length_per_dimnb.float64[3]

Cell length in each dimension of system grid

Ncellnb.int64

Total number of cells of system grid

AABB_centersnb.float64[:,3]

Center of each system grid cell

CellListCellListMesh.class_type.instance_type

Linked cell list of mesh triangles

particle_typesnb.types.DictType(key_typePT, value_typePT)

Particle types

particle_id_to_namenb.typeof(np.zeros(10, dtype = ‘U20’))

Returns the name of particle type by its id

ntypesnb.uint64

Total number of particle types

ntypes_rbnb.uint64

Total number of molecule types

kBnb.float64

Boltzmann constant

Tempnb.float64

Temperature of the system

etanb.float64

Viscosity

max_reactions_per_particlenb.int64

Maximum number of reactions that can be defined per particle (Default = 20)

reactions_leftnb.int64

Number of reactions left

molecule_typesnb.types.DictType(key_mt, value_mt)

Information about each molecule type

Reactions_Dictnb.types.DictType(keytype_RD, valuetype_RD)

Dictionary of reactions defined

reaction_idnb.int64

Current reaction id

up_reaction_idint64[:]

Reaction ids of all unimolecular particle reactions

um_reaction_idint64[:]

Reaction ids of all unimolecular molecule reactions

interaction_argsnb.typeof(np.zeros((1,1), dtype = item_t_inter))

Arguments for each partcile-particle interaction

interaction_IDsnb.types.DictType(nb.types.string, nb.int64)

Interaction ids of all partcile-particle interactions

reaction_argsnb.typeof(np.zeros((1,1), dtype = item_t_react))

Arguments for each particle-particle reaction (including fusion and enzymatic reactions)

molecule_id_to_namenb.typeof(np.zeros(10, dtype = ‘U20’))

Returns the name of molecule type by its id

mesh_scalenb.float64

Scaling factor for the compartment meshes

meshnb.types.boolean

True if mesh compartments have been defined

boundary_conditionnb.types.string (‘periodic’, ‘repulsive’, ‘fixed concentration’)

Type of boundary condition

boundary_condition_idnb.int64

Type id of boundary condition

interactionsnb.types.boolean

True if any interaction between particles has been defined

avogadronb.float64 (6.02214076e23)

Avogadros constant

kbtnb.float64

Boltzmann constant * Temperature/Avogadros constant

wall_forcenb.float64

Force constant for repulsive walls (particle-boundary and particle-triangle interactions)

barostatbarostat_lt

Berendsen barostat

seednb.int64

Random seed of simulation if defined by user

mol_type_idnb.types.DictType(nb.types.string, nb.int64)

Returns the id of molecule type by its name

react_type_idnb.types.DictType(nb.types.string, nb.int64)

Returns the id of reaction type by its name

virial_scalarnb.float64[:]

Virial (scalar)

virial_scalar_Wallnb.float64[:]

Virial (scalar) of particle-wall interaction

virial_tensornb.float64[:,:,:]

Virial tensor

virial_tensor_Wallnb.float64[:,:,:]

Virial tensor of particle-wall interaction

Pressurenb.float64[:]

Current pressure

N_bondsnb.int64[:,:]

Current number of bonds per particle type pair

namenb.types.string (‘System’)

Name of compartment

idnb.int64 (0)

ID of compartment

countnb.int64

Counter attribute (for debugging purposes)

current_stepnb.int64

Current simulation time step

time_stampnb.uint64

Gives the current number of rays that have been cast for collision detection

interaction_definednb.types.boolean[:,:]

True if any interaction or reaction between the two particle types has been defined

pair_interactionnb.types.boolean[:]

True if any pair-interaction has been defined for this particle type

length_units_prefixnb.types.DictType(nb.types.string, nb.float64)

Length unit prefix (Default = 1e6 (‘micrometer’))

time_units_prefixnb.types.DictType(nb.types.string, nb.float64)

Time unit prefix (Default = 1e0 (‘s’))

max_cells_per_dimnb.int64

Maximum number of cells of system grid in each dimension (Default = 50)

dx_randrnd.Interpolate.class_type.instance_type

Interpolation of inverse cumulative probability function describing the displacement of molecules along a planes normal after crossing by diffusion.

compartments_idnb.types.DictType(nb.types.string, nb.int64)

Returns the compartment id by its name

compartments_namenb.types.ListType(nb.types.string)

Returns the compartment name by its id

Methods

add_system_grid()	Calculates properties of the system cell grid.
add_barostat_berendsen(Tau_P, P0, start)	Adds a berendsen barostat to the system.
update_rb_barostat(self, RBs, Particles)	Updates the rigid bead molecule positions according to the scaling factor \(\mu\) of the Berendsen barostat.
register_molecule_type(molecule_name, particle_pos, particle_types, collision_type = 0)	Regsiters a molecule type.
set_diffusion_tensor(molecule_name, D_tt, D_rr, mu_tb, mu_rb, mu_tb_sqrt, mu_rb_sqrt)	Sets the diffusion tensor for a molecule type.
fixed_concentration_at_boundary(molecule_name, C, comp_name, location)	Sets the fixed concentration boundary given a concentartion for each defined molecule type and for each compartents volume and surface.
register_particle_type(Type, radius)	Registers a particle type.
set_compartments(comp_name, triangle_ids)	Sets a compartment given its triangle ids.
add_border_3d(triangle_ids_border)	Registers the border faces of an intersection between compartments and the simulation box (Needed for fixed concentration boundary condition)
add_mesh(vertices, triangles, mesh_scale = 1.0, box_triangle_ids = None, triangle_ids_transparent = None)	Set up the mesh for the mesh compartments.
add_edges()	Finds the edges of each triangle and adds them to the Mesh struct.
add_neighbours()	Finds the neighbours of each triangle and adds them to the Mesh struct.
get_AABB_centers()	Calculates the center of each cell of the system grid.
create_CellList()	Creates the linked cell list for the triangulated meshes (Used during collision detection).
..	Particle interactions
add_interaction(self, interaction_type, type1, type2, parameters, bond = False, breakable = False)
..	Uniparticle reactions
add_up_reaction(self, reaction_type, educt_type, rate, product_types, product_loc, product_direction, radius)
..	Unimolecular reactions
add_um_reaction(self, reaction_type, educt_type, rate, product_types, product_loc, product_direction, radius)
..	Biparticle reactions
add_bp_reaction(self, reaction_type, educt_types, product_types, rate, radius, interaction_type = None, interaction_parameters = None)
..	Bimolecular reactions
add_bm_reaction(self, reaction_type, educt_types, product_types, particle_pairs, pair_rates, pair_Radii)

add_barostat_berendsen(Tau_P, P0, start)

Adds a berendsen barostat to the system.

Parameters

Tau_Pfloat

Time constant of berendsen barostat

P0float

Target pressure of berendsen barostat

startint

Simulation step from which to start simulating in NPT ensemble.

Notes

The Berendsen barostat is frequently used in molecular dynamics simulations due to its simplicity and fast equilibration. While the Berendsen barostat results in the correct particle density, it does not correctly sample from the (NPT) thermodynamical ensemble, i.e. the fluctuation in the pressure are not 100% accurate. Therefore, one needs to be carefull when to use this Barostat. However, if the barostat is used in the preparation step but not during sampling, e.g. when doing LLPs simualtions where the barostat is used to relax the system’s shape but afterwards one simulates in the NVT ensemble Muller et al. [32], this should not be a problem. Note: Baruffi et al. [33] introduced a barostat for overdamped Langevin dynamics that settles to the correct thermodynamic equilibrium distribution.

add_bm_reaction(reaction_type, educt_types, product_types, particle_pairs, pair_rates, pair_Radii, placement_factor=0.5)

Registeres a Bi-Molecular (BM) reaction for an educt molecule.

Parameters

reaction_typestr {‘fusion’, ‘enzymatic_mol’}

Name of the reaction type

educt_typearray_like

Names of the educt molecule types

particle_pairsarray_like

List of educt particle pairs belonging to each of the two educt molecules

pair_ratesarray_like

List of reaction rates for each particle pair

pair_Radiiarray_like

List of reaction radii for each particle pair

placement_factorfloat64 [0,1]

Only for fusion reactions: Affects where between the educt molecules the product molecule is placed. A factor of 0.5 means in the middle. A smaller factor will shift the product towards the first educt, a larger value towards the second educt molecule. A value different from 0.5 may increase accuracy in the case of a strong size difference between the two educts (a strong difference in the diffusion tensors). Default = 0.5

Notes

Bi-molecular reactions are evaluated using the Doi scheme. Thereby, two educt particles react with a reaction rate k if the inter-particle distance is below the reaction radius. PyRID only evaluates the eucludean distances between particles but not molecule centers. Thereby, also for bi-molecular reactions educt particles need to be set. Since molecules can consist of several particles, a list of particle pairs can be passed and for each particle pair a reaction rate and radius is defined. If one of the pair particle reactions is successfull, the corresponding molecule that the particle belongs to is converted as defined by the reaction type. If the moelcule cosnsists of a single particle, the reaction will be handled just as the bi-particle reactions. The user needs to take care of calculating accurate reaction radii and rates such that the desired reaction dynamics are simulated, which can become difficult for large molecules. Therefore it is sometimes more convenient to place an interaction free particle at the origin of the molecule and execute the bimolecular reaction via this particle. This method is equivalent to the case where one would evaluate reactions based on the distance between molecule centers. As for the uni-moelcular reactions, several reaction paths can be defined. Each reaction path may involve different particle pairs of the educt molecules. On a side note: Thereby, orientation dependent bimolecular reactions are in principle possibile. However, since the reaction probability does not scale with the distance, directionality is only really accounted for if the reaction radius is small compared to the molecule size.

add_border_3d(triangle_ids_border)

Registers the border faces of an intersection between compartments and the simulation box (Needed for fixed concentration boundary condition)

Parameters

triangle_ids_borderint64[:]

IDs of the triangles making up the compartment-simulation box border.

add_bp_reaction(reaction_type, educt_types, product_types, rate, radius, interaction_type=None, interaction_parameters=None)

Registeres a Bi-Particle (BP) reaction for an educt particle.

Parameters

reaction_typestr {‘bind’, ‘absorption’, ‘enzymatic’}

Name of the reaction type

educt_typearray_like

Names of the educt particle types

product_typesarray_like

List of products (only used for binding reactions)

ratefloat64

Reaction rate

radiusfloat64

Reaction radius.

interaction_typestr

Name of the energy potential function in case ‘reaction_type’ is a binding reaction.

interaction_parametersdict

Parameters for the interaction function in case ‘reaction_type’ is a binding reaction.

Notes

Bi-particle reactions are evaluated using the Doi scheme. Thereby, two educt particles react with a reaction rate k if the inter-particle distance is below the reaction radius.

Binding Reactions

A particle pair will only enter a bound state if the binding reaction was succesfull. Otherwise two particles won’t experience any pairwise interaction force! A reaction occurs if two particles for whose type a binding reaction has been defined are within the reaction radius. The reaction probability is evaluated as:

\[1-exp(\lambda \cdot \Delta t),\]

where \(\lambda\) is the reaction rate and \(\Delta t\) is the integration time step.

add_edges()

Finds the edges of each triangle and adds them to the Mesh struct.

Notes

Edges are sorted in counter clockwise order:

edge: [0,1]

[1,2]

[0,2]

add_interaction(interaction_type, type1, type2, parameters, bond=False, breakable=False)

Assigns a particle-particle interaction potential to a particle type pair.

Parameters

interaction_typestr {‘harmonic_repulsion’, ‘PHS’, ‘harmonic_attraction’, ‘CSW’}

Type of interaction potential

type1str

Names of the educt particle type 1

type2str

Names of the educt particle type 2

parametersdict

Parameters for the interaction potential function.

bondbool

True if interaction is initialized as part of a bindign reaction. (Default: False)

breakablebool

If set to True, bonds that have been formed by a binding reaction are deleted if the interparticle distance is above the cutoff radius of the interaction potenial.

Notes

Bi-particle reactions are evaluated using the Doi scheme. Thereby, two educt particles react with a reaction rate k if the inter-particle distance is below the reaction radius.

add_mesh(vertices, triangles, mesh_scale=1.0, box_triangle_ids=None, triangle_ids_transparent=None, adapt_box=False)

Set up the mesh for the mesh compartments.

Parameters

verticesfloat64[:,3]

Position of mesh vertices in each dimension.

trianglesint64[:,3]

Indices of vertices that make up each triangle (pointers to vertices)

mesh_scalefloat

Scaling factor to scale up or down the mesh (Default = 1.0).

box_triangle_idsint64[:]

Triangle ids of the faces that make up the simulation box border (To be excluded from self.triangle_ids)

triangle_ids_transparentint64[:]

Triangle ids of any transparent faces (To be excluded from self.triangle_ids)

adapt_boxboolean

If true, the simulation box size is adapted such that the mesh fits exactly in the simulation box. Default = False

Notes

A Compartment is defined by a triangulated manifold mesh (a mesh without holes and disconnected vertices or edges), i.e. it has no gaps and separates the space on the inside of the surface from the space outside. It also looks like a surface everywhere on the mesh Shirley et al. [34]. The mesh is stored as a shared vertex mesh. This data structure has triangles which point to vertices which contain the vertex data:

triangles: Three vectors per triangle

0 | (0,1,2)

1 | (1,3,2)

2 | (0,3,1)

vertices: One vector per vertex

0 | (a_x,a_y,a_z)

1 | (b_x,b_y,b_z)

2 | (c_x,c_y,c_z)

3 | (d_x,d_y,d_z)

During raymarching on a mesh surface, we need to find the neighbor triangle given some edge i that has been crossed. We may save, for each triangle, the corresponding edges in an array:

0 | (0,1), (0,2), (1,2)

1 | …

In addition, we save, for each edge of a triangle the corresponding neighbor in another array:

0 | (1, 4, 2)

1 | …

Also, we need to rotate the particles between local coordinate systems of neighbouring triangles, as particles cross an edge. We may precalculate the necessary parameters and save them in a datastructure similar to the above:

0 | (params_0->1, params_0->4, params_0->2)

1 | …

params_i->j : [sinPhi, cosPhi, a_n]

We may also want to prcompute the rotation parameters for each triangle pait to speed up the coordiante transformation when a particle hops between two triangles. We may take the same data structure as before:

0 | (rot_params[0->1], (rot_params[0->4], (rot_params[0->2])

1 | …

rot_params: array[sin, cos, an]

add_neighbours()

Finds the neighbours of each triangle and adds them to the Mesh struct.

add_system_grid()

Calculates properties of the system cell grid.

add_um_reaction(reaction_type, educt_type, rate, product_types, product_loc, product_direction, radius)

Registeres a Uni-Molecular (UM) reaction for an educt Molecule.

Parameters

reaction_typestr {‘production’, ‘fission’, ‘conversion_mol’, ‘decay’}

Name of the reaction type

educt_typestr

Name of the educt molecule type

ratefloat64

Reaction rate

product_typesarray_like

List of products

product_locarray_like

List specifying whether the individual products are volume (0) or surface (1) molecules

product_directionarray_like

Specifies whether, in case of a surface molecule releasing products into the volume, the realease occurs inside (-1) or outside (1) the mesh compartment.

radiusfloat64

The radius within which molecules are distributed upon a successfull release reaction.

Notes

Uni-molecular reactions are evaluated using the Gillespie stochastic simulation algorithm. Thereby, for a molecule type, the time point of the next uni-molecular reaction occuring is calculated in advance and executed when the simualtion reaches the respective time point. A uni-molecular reaction can consist of several reaction paths. Each molecule type has exactly one uni-molecular reaction id registered for which the time point is drawn randomly based on the total reaction rate which is the sum of all path reaction rates. At the time point of the reaction occuring, a path is randomly drawn from the set of possible paths weighted by the individual path reaction rates. To illustrate the above by an example: For Molecule ‘A’ a we may define two convcersion reactions where A -k1-> B, A -k2-> C. In addition, we define a fission reaction A -k3-> D+E. The next reaction will occur at a time point T that is calculated based on the total rate k1+k2+k3. At time T, one of the three reactions is executed. The probability for each reaction to be selected being k1/(k1+k2+k3), k2/(k1+k2+k3), k3/(k1+k2+k3).

add_up_reaction(reaction_type, educt_type, rate, product_types, product_loc, product_direction, radius)

Registeres a Uni-Particle (UP) reaction for an educt particle.

Parameters

reaction_typestr {‘relase’, ‘conversion’}

Name of the reaction type

educt_typestr

Name of the educt particle type

ratefloat64

Reaction rate

product_typesarray_like

List of products

product_locarray_like

List specifying whether the individual products are volume (0) or surface (1) molecules

product_directionarray_like

Specifies whether, in case of a surface particle releasing products into the volume, the realease occurs inside (-1) or outside (1) the mesh compartment.

radiusfloat64

The radius within which molecules are distributed upon a successfull release reaction.

Notes

Uni-particle reactions are evaluated using a variant of the Gillespie stochastic simulation algorithm. Thereby, for a particle type, the time point of the next uni-particle reaction occuring is calculated in advance and executed when the simualtion reaches the respective time point. For unimolecular reactions we can draw the time point of the next reaction from a distribution using a variant of the Gillespie Stochastic Simulation Algorithm (SSA) Stiles and Bartol [11], Erban et al. [12]. For a single molecule having \(n\) possible transition reactions/reaction paths each having a reaction rate \(k_i\), let \(k_t = \sum_i^n k_i\) be the total reaction rate.

Now, let \(\rho(\tau) d\tau\) be the probability that the next reaction occurs within \([t+\tau,t+\tau+d\tau)\) and let \(g(\tau)\) be the probability that no reaction occurs within \([t,t+\tau)\). The probability that a reaction occurs within the time interval \(d\tau\) is simply \(k_t d\tau\). Thereby

\[\rho(\tau) d\tau = g(\tau) k_t d\tau\]

where \(g(\tau) = e^{-k_t \tau}\) (see also Poisson process). From the above equation we easily find \(P(\tau) = 1-e^{-k_t \tau}\) by integration. To sample from this distribution, we can use the inverse distribution function.

\[\tau = P^{-1}(U)\]

where \(U\) i uniformly distributed in \((0,1)\). Given that \(U = P(\tau) = 1-e^{-k_t \tau}\), we find \(P^{-1}(U) = \frac{-log(1-U)}{k_t}\). Since \(U\) is uniformly distributed in 0,1, so is \(1-U\). Thereby, we can draw the time point of the next reaction from:

\[\tau = \frac{1}{k_t} \ln\Big[\frac{1}{U}\Big],\]

With the above method, we accurately sample from the distribution of expected molecule lifetimes \(\rho(\tau) = k_t e^{-k_t \tau}\).

A uni-particle reaction can consist of several reaction paths. Each particle type has exactly one uni-particle reaction id registered for which the time point is drawn randomly based on the total reaction rate which is the sum of all path reaction rates. At the time point of the reaction occuring, a path is randomly drawn from the set of possible paths weighted by the individual path reaction rates. To get the reaction path we compare a second random number, uninformly distributed in (0,k_t), with the cumulative set of reaction rates \((k_1, k_1+k_2, ..., k_t)\). The comparison is done via a bisection algorithm (see rand_util.random_util.bisect_right()).

To illustrate the above by an example: For particle ‘A’ a we may define two convcersion reactions where A -k1-> B, A -k2-> C. In addition, we define a release reaction A -k3-> 2*D+5*E. The next reaction will occur at a time point T that is calculated based on the total rate k1+k2+k3. At time T, one of the three reactions is executed. The probability for each reaction to be selected being k1/(k1+k2+k3), k2/(k1+k2+k3), k3/(k1+k2+k3).

create_cell_list()

Creates the linked cell list for the triangulated meshes by calling system_util.CellList_util.create_cell_list_mesh() (Used during collision detection).

fixed_concentration_at_boundary(molecule_name, C, comp_name, location)

Sets the fixed concentration boundary given a concentartion for each defined molecule type and for each compartents volume and surface.

Parameters

molecule_namestr

Name of the molecule

Cfloat

concentration of molecule type outside simulation box (‘virtual molecules’)

comp_namestr

Name of the compartment

locationint {0,1}

Location of the molecule. 0: Volume, 1: Surface

Notes

fixed_concentration_at_boundary() calculates some properties that are necessary to properly distribute molecules inside the simulation box that hit the simulation box boundary from the outside (Thereby, ‘virtual molecules’ become real molecules in our simualtion). The number of hits per time step a boundary of area A experiences is \(N = l_{perp}*A*C\). Where \(C\) is the concentration in molecules per volume and \(l_{perp}\) is the average net displacement in one tiem step towards or away from any plane, where \(l_{perp} = \sqrt{(4*D*\Delta t/\pi)}\) Stiles and Bartol [11].

get_AABB_centers()

Calculates the center of each cell of the system grid.

register_molecule_type(molecule_name, particle_pos, particle_types, collision_type=0, h_membrane=None)

Regsiters a molecule type.

Parameters

molecule_namestr

Name of the molecule.

particle_posint[:,3]

Positions of each of the molecule’s particles.

particle_typesarray_like (dtype = ‘U20’)

Array of names for each of the molecule’s particles.

collision_typeint {0,1}

Collision type of the molecule (Default = 0). 0: Collisions with compartments are handled via interaction force. 1: Collisions with compartments are handled by raytracing algorithm (recommended for molcules whose diffusion speed and/or interaction radius is small compared to the chosen integration time step).

h_membranefloat64

Distance of the molecules´ center of diffusion from the membrane surface. Default = None

register_particle_type(Type, radius)

Registers a particle type.

Parameters

Typestr

Name of the particle type.

radiusfloat

radius of the particle.

set_compartments(comp_name, triangle_ids)

Sets a compartment given its triangle ids.

Parameters

comp_namestr

Name of the compartment.

triangle_idsint64[:]

IDs of the triangles making up the compartment mesh.

set_diffusion_tensor(molecule_name, D_tt, D_rr, mu_tb, mu_rb, mu_tb_sqrt, mu_rb_sqrt)

Sets the diffusion tensor for a given molecule type.

Parameters

molecule_namestr

Some Information

D_ttfloat[3,3]

Translational diffusion tensor ()

D_rrfloat[3,3]

Rotational diffusion tensor

mu_tbfloat[3,3]

Translational mobility tensor

mu_rbfloat[3,3]

Rotational mobility tensor

mu_tb_sqrtfloat[3,3]

Square root of translational mobility tensor

mu_rb_sqrtfloat[3,3]

Square root of rotational mobility tensor

Notes

PyRID offers a module for the calculation of the diffusion tensors of rigid bead models (see molecules_util.hydro_util). It uses the Kirkwood-Riseman calculation with modified Oseen tensor as introduced by Carrasco and de la Torre [2], Carrasco and de la Torre [16]. The same method has been implemented in the HYDRO++ tool, which is free to use but not open source! As such, you may also use HYDRO++ for the calculation of the diffusion tensors. Or even better, use both and tell me if you you get inconsistent results ;) !

Note

Diffusion tensor should be a real positive semidefinite matrix!