diff --git a/simularium_readdy_models/actin/actin_generator.py b/simularium_readdy_models/actin/actin_generator.py index 9ffa3c7..ea3fd28 100644 --- a/simularium_readdy_models/actin/actin_generator.py +++ b/simularium_readdy_models/actin/actin_generator.py @@ -449,6 +449,7 @@ def _get_main_monomers_for_fiber( pointed_actin_number, particles={}, longitudinal_bonds=True, + barbed_binding_site=False, ): """ get the main actins for a fiber (i.e. no branches or arps). @@ -549,8 +550,14 @@ def _get_main_monomers_for_fiber( # remove "mid" from second actin ActinGenerator._remove_mid_from_actin(particle_ids[1], particles) actin_arp_ids = None + if barbed_binding_site and len(particle_ids) > 1: + particles = ActinGenerator._set_particle_type_name( + "actin#barbed_ATP_" , + particle_ids[len(particle_ids) - 2], + particles, + ) particles = ActinGenerator._set_particle_type_name( - "actin#barbed_ATP_", + "actin#barbed_ATP_" if not barbed_binding_site else "binding_site#" , particle_ids[len(particle_ids) - 1], particles, ) @@ -565,6 +572,7 @@ def _get_monomers_for_fiber( pointed_actin_number, particles={}, longitudinal_bonds=True, + barbed_binding_site=False, ): """ get the main actins for a fiber as well as any bound arps and daughter fibers. @@ -581,6 +589,7 @@ def _get_monomers_for_fiber( pointed_actin_number, particles, longitudinal_bonds, + barbed_binding_site, ) daughter_particle_ids = [] all_actin_arp_ids = [] @@ -800,6 +809,7 @@ def get_monomers( use_uuids=True, start_normal=None, longitudinal_bonds=True, + barbed_binding_site=False, ): """ get all the monomer data for the (branched) fibers in fibers_data. @@ -836,6 +846,7 @@ def get_monomers( np.zeros(3), 1, longitudinal_bonds=longitudinal_bonds, + barbed_binding_site=barbed_binding_site, ) result["topologies"][ActinGenerator._get_next_monomer_id()] = { "type_name": "Actin-Polymer", diff --git a/simularium_readdy_models/actin/actin_simulation.py b/simularium_readdy_models/actin/actin_simulation.py index a0826c7..042c4fa 100644 --- a/simularium_readdy_models/actin/actin_simulation.py +++ b/simularium_readdy_models/actin/actin_simulation.py @@ -4,7 +4,13 @@ import numpy as np import readdy -from ..common import ReaddyUtil +from ..common import ( + ReaddyUtil, + add_membrane_particle_types, + add_membrane_constraints, + init_membrane, + all_membrane_particle_types +) from .actin_structure import ActinStructure from .actin_util import ActinUtil @@ -110,7 +116,14 @@ def add_particle_types(self): self.actin_util.add_actin_types(self.system, actin_diffCoeff) self.actin_util.add_arp23_types(self.system, arp23_diffCoeff) self.actin_util.add_cap_types(self.system, cap_diffCoeff) - self.system.add_species("obstacle", 0.0) + self.system.topologies.add_type("Obstacle") + self.system.add_topology_species("obstacle", self._parameter("obstacle_diff_coeff")) + if self._parameter("add_membrane"): + add_membrane_particle_types( + self.system, self._parameter("membrane_particle_radius"), temperature, viscosity + ) + if self._parameter("barbed_binding_site"): + self.actin_util.add_binding_site_types(self.system, actin_diffCoeff) def add_constraints(self): """ @@ -120,6 +133,7 @@ def add_constraints(self): util = ReaddyUtil() longitudinal_bonds = bool(self._parameter("longitudinal_bonds")) only_linear_actin = bool(self._parameter("only_linear_actin_constraints")) + actin_constraints = bool(self._parameter("actin_constraints")) # force constants actin_angle_force_constant = float(self._parameter("angles_force_constant")) actin_dihedral_force_constant = float( @@ -132,16 +146,17 @@ def add_constraints(self): longitudinal_bonds, float(self._parameter("bonds_force_multiplier")), ) - self.actin_util.add_filament_twist_angles( - actin_angle_force_constant, self.system, util, longitudinal_bonds - ) - self.actin_util.add_filament_twist_dihedrals( - actin_dihedral_force_constant, - self.system, - util, - longitudinal_bonds, - only_linear_actin, - ) + if actin_constraints: + self.actin_util.add_filament_twist_angles( + actin_angle_force_constant, self.system, util, longitudinal_bonds + ) + self.actin_util.add_filament_twist_dihedrals( + actin_dihedral_force_constant, + self.system, + util, + longitudinal_bonds, + only_linear_actin, + ) if not only_linear_actin: # branch junction self.actin_util.add_branch_bonds(self.system, util) @@ -170,8 +185,41 @@ def add_constraints(self): actin_actin_repulsion_potentials=True, longitudinal_bonds=longitudinal_bonds, ) + self.actin_util.add_repulsions_with_actin( + ["obstacle"], + self._parameter("obstacle_radius"), + ActinUtil.DEFAULT_FORCE_CONSTANT, + self.system, + util + ) # box potentials self.actin_util.add_monomer_box_potentials(self.system) + self.actin_util.add_obstacle_box_potential(self.system) + self.actin_util.add_extra_box(self.system) + # membrane + if self._parameter("add_membrane"): + add_membrane_constraints( + self.system, + np.array([ + float(self._parameter("membrane_center_x")), + float(self._parameter("membrane_center_y")), + float(self._parameter("membrane_center_z")), + ]), + np.array([ + float(self._parameter("membrane_size_x")), + float(self._parameter("membrane_size_y")), + float(self._parameter("membrane_size_z")), + ]), + self._parameter("membrane_particle_radius"), + self._parameter("box_size") + ) + self.actin_util.add_repulsions_with_actin( + all_membrane_particle_types(), + self._parameter("membrane_particle_radius"), + ActinUtil.DEFAULT_FORCE_CONSTANT, + self.system, + util + ) def add_reactions(self): """ @@ -198,6 +246,8 @@ def add_reactions(self): self.actin_util.add_cap_unbind_reaction(self.system) if self.do_pointed_end_translation(): self.actin_util.add_translate_reaction(self.system) + if self._parameter("position_obstacle_stride") > 0: + self.actin_util.add_position_obstacle_reaction(self.system) def do_pointed_end_translation(self): result = self._parameter("displace_pointed_end_tangent") or self._parameter( @@ -330,19 +380,44 @@ def add_obstacles(self): """ Add obstacle particles. """ + if not self._parameter("add_obstacles"): + return n = 0 while f"obstacle{n}_position_x" in self.parameters: - self.simulation.add_particle( - type="obstacle", - position=[ + self.simulation.add_topology( + "Obstacle", + ["obstacle"], + np.array([[ float(self._parameter(f"obstacle{n}_position_x")), float(self._parameter(f"obstacle{n}_position_y")), float(self._parameter(f"obstacle{n}_position_z")), - ], + ]]) ) n += 1 if n > 0: print(f"Added {n} obstacle(s).") + + def add_membrane(self): + """ + Add membrane types and particles. + """ + if not self._parameter("add_membrane"): + return + init_membrane( + self.simulation, + np.array([ + float(self._parameter("membrane_center_x")), + float(self._parameter("membrane_center_y")), + float(self._parameter("membrane_center_z")), + ]), + np.array([ + float(self._parameter("membrane_size_x")), + float(self._parameter("membrane_size_y")), + float(self._parameter("membrane_size_z")), + ]), + self._parameter("membrane_particle_radius"), + self._parameter("box_size") + ) def add_crystal_structure_monomers(self): """ diff --git a/simularium_readdy_models/actin/actin_util.py b/simularium_readdy_models/actin/actin_util.py index a8b3824..c354aad 100644 --- a/simularium_readdy_models/actin/actin_util.py +++ b/simularium_readdy_models/actin/actin_util.py @@ -33,6 +33,10 @@ def set_displacements(d): pointed_monomer_positions = [] +obstacle_time_index = 0 +obstacle_position = np.array([41., 0, 0]) # TODO set from another simulator etc + + class ActinUtil: DEFAULT_FORCE_CONSTANT = 250.0 @@ -98,6 +102,9 @@ def __init__(self, parameters, displacements=None): "use_box_arp": False, "use_box_cap": False, "obstacle_radius": 35.0, + "obstacle_diff_coeff": 0.0, + "use_box_obstacle": False, + "position_obstacle_stride": 0, "n_fixed_monomers_pointed": 0, "n_fixed_monomers_barbed": 0, "displace_pointed_end_tangent": False, @@ -117,6 +124,13 @@ def __init__(self, parameters, displacements=None): "bonds_force_multiplier": 0.2, "angles_force_constant": 1000.0, "dihedrals_force_constant": 1000.0, + "add_membrane": False, + "add_obstacles": False, + "actin_constraints": True, + "add_monomer_box_potentials": False, + "add_extra_box": False, + "barbed_binding_site": False, + "binding_site_reaction_distance": 0.1, } @staticmethod @@ -134,6 +148,21 @@ def get_new_vertex(topology): ) return results[0] + @staticmethod + def get_binding_site_vertex(topology): + """ + Get the barbed end binding site vertex. + """ + results = ReaddyUtil.get_vertices_of_type( + topology, "binding_site", exact_match=False, error_msg="Failed to find binding site vertex" + ) + if len(results) > 1: + raise Exception( + f"Found more than one binding site vertex\n" + f"{ReaddyUtil.topology_to_string(topology)}" + ) + return results[0] + @staticmethod def get_new_arp23(topology): """ @@ -164,11 +193,6 @@ def get_actin_number(topology, vertex, offset): (i.e. return 3 for type = "actin#ATP_1" and offset = -1). """ pt = topology.particle_type_of_vertex(vertex) - if "actin" not in pt: - raise Exception( - f"Failed to get actin number: {pt} is not actin\n" - f"{ReaddyUtil.topology_to_string(topology)}" - ) return ReaddyUtil.calculate_polymer_number( int(pt[-1]), offset, ActinUtil.n_polymer_numbers() ) @@ -252,14 +276,15 @@ def get_next_actin(topology, v_actin, direction, error_if_not_found=False): get the next actin toward the pointed or barbed direction. """ n = ActinUtil.get_actin_number(topology, v_actin, direction) - end_type = "barbed" if direction > 0 else "pointed" actin_types = [ f"actin#ATP_{n}", f"actin#{n}", f"actin#mid_ATP_{n}", f"actin#mid_{n}", - f"actin#{end_type}_ATP_{n}", - f"actin#{end_type}_{n}", + f"actin#barbed_ATP_{n}", + f"actin#barbed_{n}", + f"actin#pointed_ATP_{n}", + f"actin#pointed_{n}", ] if direction < 0 and n == 1: actin_types += ["actin#branch_1", "actin#branch_ATP_1"] @@ -385,8 +410,9 @@ def set_end_vertex_position(topology, recipe, v_new, barbed): else ActinStructure.mother1_to_mother_vector() ) at_branch = False + direction = -1 if barbed else 1 vertices.append( - ReaddyUtil.get_neighbor_of_type(topology, v_new, "actin", False) + ActinUtil.get_next_actin(topology, v_new, direction, False) ) if vertices[0] is None: ( @@ -396,9 +422,7 @@ def set_end_vertex_position(topology, recipe, v_new, barbed): at_branch = True else: vertices.append( - ReaddyUtil.get_neighbor_of_type( - topology, vertices[0], "actin", False, [v_new] - ) + ActinUtil.get_next_actin(topology, vertices[0], direction, False) ) if vertices[1] is None: ( @@ -410,9 +434,7 @@ def set_end_vertex_position(topology, recipe, v_new, barbed): at_branch = True else: vertices.append( - ReaddyUtil.get_neighbor_of_type( - topology, vertices[1], "actin", False, [vertices[0]] - ) + ActinUtil.get_next_actin(topology, vertices[1], direction, False) ) if vertices[2] is None: ( @@ -1013,8 +1035,33 @@ def reaction_function_finish_barbed_grow(topology): """ reaction function for the barbed end growing. """ + if parameters["barbed_binding_site"]: + return ActinUtil.do_finish_binding_site_grow(topology) return ActinUtil.do_finish_grow(topology, True) + @staticmethod + def do_finish_binding_site_grow(topology): + """ + reaction function for the barbed end growing with a binding site. + """ + recipe = readdy.StructuralReactionRecipe(topology) + end_type = "barbed" + if parameters["verbose"]: + print("Grow " + end_type) + v_bs = ActinUtil.get_binding_site_vertex(topology) + v_barbed = ReaddyUtil.get_first_neighbor( + topology, + v_bs, + [], + error_msg=f"Failed to find neighbor of new {end_type} end", + ) + v_neighbor = ActinUtil.get_next_actin(topology, v_barbed, -1, True) + ActinUtil.set_end_vertex_position(topology, recipe, v_bs, True) + recipe.add_edge(v_bs, v_neighbor) + ReaddyUtil.set_flags(topology, recipe, v_neighbor, ["mid"], ["barbed"], True) + recipe.change_topology_type("Actin-Polymer") + return recipe + @staticmethod def reaction_function_finish_arp_bind(topology): """ @@ -1483,6 +1530,23 @@ def reaction_function_translate(topology): time_index += 1 return recipe + @staticmethod + def reaction_function_position_obstacle(topology): + """ + reaction function to set position of obstacle particle. + """ + global obstacle_time_index + recipe = readdy.StructuralReactionRecipe(topology) + if obstacle_time_index > 0 and obstacle_time_index % parameters["position_obstacle_stride"] != 0: + obstacle_time_index += 1 + return recipe + if parameters["verbose"]: + print("Translate obstacle") + v = topology.graph.get_vertices()[0] # there should only be one particle in topology + recipe.change_particle_position(v, obstacle_position) + obstacle_time_index += 1 + return recipe + @staticmethod def get_all_actin_particle_types(): """ @@ -1585,6 +1649,16 @@ def get_all_cap_particle_types(): "cap#bound", ] + @staticmethod + def get_all_binding_site_particle_types(): + """ + get particle types for barbed binding sites. + """ + result = [] + for i in ActinUtil.polymer_number_range(): + result += [f"binding_site#{i}"] + return result + @staticmethod def get_all_particle_types(): """ @@ -1656,6 +1730,15 @@ def add_cap_types(system, diffCoeff): ActinUtil.get_all_cap_particle_types(), system, diffCoeff ) + @staticmethod + def add_binding_site_types(system, diffCoeff): + """ + add particle types for binding sites (only implemented for barbed end, ATP). + """ + ActinUtil.add_particle_types( + ActinUtil.get_all_binding_site_particle_types(), system, diffCoeff + ) + @staticmethod def add_bonds_between_actins(system, util, longitudinal_bonds, force_multiplier): """ @@ -1666,82 +1749,90 @@ def add_bonds_between_actins(system, util, longitudinal_bonds, force_multiplier) n_polymer_numbers = ActinUtil.n_polymer_numbers() lat_force_constant = force_multiplier * 968.2 # kJ / mol / nm^2 long_force_constant = force_multiplier * 1437.5 # kJ / mol / nm^2 + pointed_actins = [ + "actin#", + "actin#ATP_", + "actin#mid_", + "actin#mid_ATP_", + "actin#pointed_", + "actin#pointed_ATP_", + "actin#fixed_", + "actin#fixed_ATP_", + "actin#mid_fixed_", + "actin#mid_fixed_ATP_", + "actin#pointed_fixed_", + "actin#pointed_fixed_ATP_", + ] + barbed_actins = [ + "actin#", + "actin#ATP_", + "actin#mid_", + "actin#mid_ATP_", + "actin#barbed_", + "actin#barbed_ATP_", + "actin#fixed_", + "actin#fixed_ATP_", + "actin#mid_fixed_", + "actin#mid_fixed_ATP_", + "actin#fixed_barbed_", + "actin#fixed_barbed_ATP_", + ] # lateral actin-actin bond util.add_polymer_bond_1D( - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#pointed_", - "actin#pointed_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#pointed_fixed_", - "actin#pointed_fixed_ATP_", - ], + pointed_actins, 0, - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#barbed_", - "actin#barbed_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#fixed_barbed_", - "actin#fixed_barbed_ATP_", - ], + barbed_actins, 1, lat_force_constant, bond_length_lat, system, n_polymer_numbers, ) + if parameters["barbed_binding_site"]: + util.add_polymer_bond_1D( + ["actin#barbed_", "actin#barbed_ATP_"], + 0, + ["binding_site#"], + 1, + lat_force_constant, + bond_length_lat, + system, + n_polymer_numbers, + ) + util.add_polymer_bond_1D( + ["actin#barbed_", "actin#barbed_ATP_"], + 0, + ["actin#barbed_", "actin#barbed_ATP_"], + 1, + lat_force_constant, + bond_length_lat, + system, + n_polymer_numbers, + ) print(f"Added lat bonds with fc = {lat_force_constant}") if longitudinal_bonds: print("Adding longitudinal bonds...") util.add_polymer_bond_1D( - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#pointed_", - "actin#pointed_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#pointed_fixed_", - "actin#pointed_fixed_ATP_", - ], + pointed_actins, 0, - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#barbed_", - "actin#barbed_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#fixed_barbed_", - "actin#fixed_barbed_ATP_", - ], + barbed_actins, 2, long_force_constant, bond_length_long, system, n_polymer_numbers, ) + if parameters["barbed_binding_site"]: + util.add_polymer_bond_1D( + pointed_actins, + 0, + ["binding_site#"], + 2, + long_force_constant, + bond_length_long, + system, + n_polymer_numbers, + ) print(f"Added long bonds with fc = {long_force_constant}") # branch actin-actin bond util.add_bond( @@ -1768,7 +1859,7 @@ def add_bonds_between_actins(system, util, longitudinal_bonds, force_multiplier) system, ) # temporary bonds - util.add_polymer_bond_1D( # temporary during growth reactions + util.add_polymer_bond_1D( [ "actin#", "actin#ATP_", @@ -1798,7 +1889,7 @@ def add_bonds_between_actins(system, util, longitudinal_bonds, force_multiplier) system, n_polymer_numbers, ) - util.add_bond( # temporary during growth reactions + util.add_bond( [ "actin#branch_1", "actin#branch_ATP_1", @@ -1820,49 +1911,52 @@ def add_filament_twist_angles(force_constant, system, util, longitudinal_bonds): """ add angles for filament twist and cohesiveness. """ + pointed_actins = [ + "actin#", + "actin#ATP_", + "actin#mid_", + "actin#mid_ATP_", + "actin#pointed_", + "actin#pointed_ATP_", + "actin#fixed_", + "actin#fixed_ATP_", + "actin#mid_fixed_", + "actin#mid_fixed_ATP_", + "actin#pointed_fixed_", + "actin#pointed_fixed_ATP_", + ] + mid_actins = [ + "actin#", + "actin#ATP_", + "actin#mid_", + "actin#mid_ATP_", + "actin#fixed_", + "actin#fixed_ATP_", + "actin#mid_fixed_", + "actin#mid_fixed_ATP_", + ] + barbed_actins = [ + "actin#", + "actin#ATP_", + "actin#mid_", + "actin#mid_ATP_", + "actin#barbed_", + "actin#barbed_ATP_", + "actin#fixed_", + "actin#fixed_ATP_", + "actin#mid_fixed_", + "actin#mid_fixed_ATP_", + "actin#fixed_barbed_", + "actin#fixed_barbed_ATP_", + ] # Lateral bond to lateral bond angle angle = ActinStructure.actin_to_actin_angle() util.add_polymer_angle_1D( - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#pointed_", - "actin#pointed_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#pointed_fixed_", - "actin#pointed_fixed_ATP_", - ], + pointed_actins, -1, - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - ], + mid_actins, 0, - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#barbed_", - "actin#barbed_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#fixed_barbed_", - "actin#fixed_barbed_ATP_", - ], + barbed_actins, 1, force_constant, angle, @@ -1898,52 +1992,30 @@ def add_filament_twist_angles(force_constant, system, util, longitudinal_bonds): angle, system, ) + if parameters["barbed_binding_site"]: + util.add_polymer_angle_1D( + pointed_actins, + -1, + ["actin#barbed_", "actin#barbed_ATP_"], + 0, + ["binding_site#"], + 1, + force_constant, + angle, + system, + ActinUtil.n_polymer_numbers(), + ) if not longitudinal_bonds: print(f"Added angles with fc = {force_constant}") return # Lateral bond to longitudinal bond angle angle = ActinStructure.actin_to_actin_angle(True, False) util.add_polymer_angle_1D( - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#pointed_", - "actin#pointed_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#pointed_fixed_", - "actin#pointed_fixed_ATP_", - ], + pointed_actins, -1, - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - ], + mid_actins, 0, - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#barbed_", - "actin#barbed_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#fixed_barbed_", - "actin#fixed_barbed_ATP_", - ], + barbed_actins, 2, force_constant, angle, @@ -1979,49 +2051,27 @@ def add_filament_twist_angles(force_constant, system, util, longitudinal_bonds): angle, system, ) + if parameters["barbed_binding_site"]: + util.add_polymer_angle_1D( + pointed_actins, + -1, + mid_actins, + 0, + ["binding_site#"], + 2, + force_constant, + angle, + system, + ActinUtil.n_polymer_numbers(), + ) # Longitudinal bond to longitudinal bond angle angle = ActinStructure.actin_to_actin_angle(False, False) util.add_polymer_angle_1D( - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#pointed_", - "actin#pointed_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#pointed_fixed_", - "actin#pointed_fixed_ATP_", - ], + pointed_actins, -2, - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - ], + mid_actins, 0, - [ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#barbed_", - "actin#barbed_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#fixed_barbed_", - "actin#fixed_barbed_ATP_", - ], + barbed_actins, 2, force_constant, angle, @@ -2057,6 +2107,19 @@ def add_filament_twist_angles(force_constant, system, util, longitudinal_bonds): angle, system, ) + if parameters["barbed_binding_site"]: + util.add_polymer_angle_1D( + pointed_actins, + -2, + mid_actins, + 0, + ["binding_site#"], + 2, + force_constant, + angle, + system, + ActinUtil.n_polymer_numbers(), + ) print(f"Added angles (incl longitudinal) with fc = {force_constant}") @staticmethod @@ -2904,6 +2967,45 @@ def add_actin_actin_repulsions( system, n_polymer_numbers, ) + util.add_polymer_repulsions_1D( + [ + "actin#free", + "actin#free_ATP", + ], + None, + [ + "actin#", + "actin#ATP_", + "actin#mid_", + "actin#mid_ATP_", + "actin#barbed_", + "actin#barbed_ATP_", + "actin#fixed_", + "actin#fixed_ATP_", + "actin#mid_fixed_", + "actin#mid_fixed_ATP_", + "actin#fixed_barbed_", + "actin#fixed_barbed_ATP_", + ], + 0, + force_constant, + ActinStructure.actin_to_actin_repulsion_distance(True), + system, + n_polymer_numbers, + ) + util.add_repulsion( + [ + "actin#free", + "actin#free_ATP", + ], + [ + "actin#free", + "actin#free_ATP", + ], + force_constant, + ActinStructure.actin_to_actin_repulsion_distance(True), + system, + ) if longitudinal_bonds: util.add_polymer_repulsions_1D( [ @@ -2968,13 +3070,6 @@ def add_repulsions( ActinUtil.add_actin_actin_repulsions( force_constant, system, util, longitudinal_bonds ) - util.add_repulsion( - actin_types, - ["obstacle"], - force_constant, - actin_radius + obstacle_radius, - system, - ) # arp2/3 util.add_repulsion( arp_types, @@ -2990,13 +3085,6 @@ def add_repulsions( arp23_radius + actin_radius, system, ) - util.add_repulsion( - arp_types, - ["obstacle"], - force_constant, - arp23_radius + obstacle_radius, - system, - ) # capping protein util.add_repulsion( cap_types, @@ -3019,11 +3107,22 @@ def add_repulsions( cap_radius + arp23_radius, system, ) + + @staticmethod + def add_repulsions_with_actin(other_types, other_radius, force_constant, system, util): + """ + Add repulsions between actin etc types and a given list of types. + """ + actin_types = ( + ActinUtil.get_all_actin_particle_types() + + ActinUtil.get_all_fixed_actin_particle_types() + ) + actin_radius = 0.5 * ActinStructure.actin_to_actin_repulsion_distance(True) util.add_repulsion( - cap_types, - ["obstacle"], + actin_types, + other_types, force_constant, - cap_radius + obstacle_radius, + actin_radius + other_radius, system, ) @@ -3062,8 +3161,10 @@ def check_add_global_box_potential(system): @staticmethod def add_monomer_box_potentials(system): """ - Confine free monomers to boxes centered at origin with extent. + Confine free monomers to boxes centered at center with extent. """ + if not bool(parameters["add_monomer_box_potentials"]): + return particle_types = { "actin": ["actin#free", "actin#free_ATP"], "arp": ["arp2#free"], @@ -3095,6 +3196,64 @@ def add_monomer_box_potentials(system): system, ) + @staticmethod + def add_obstacle_box_potential(system): + """ + Confine obstacle to a box centered at center with extent. + """ + if not parameters[f"use_box_obstacle"]: + return + center = np.array( + [ + parameters[f"obstacle_box_center_x"], + parameters[f"obstacle_box_center_y"], + parameters[f"obstacle_box_center_z"], + ] + ) + size = np.array( + [ + parameters[f"obstacle_box_size_x"], + parameters[f"obstacle_box_size_y"], + parameters[f"obstacle_box_size_z"], + ] + ) + ActinUtil.add_box_potential( + ["obstacle"], + center - 0.5 * size, + size, + ActinUtil.DEFAULT_FORCE_CONSTANT, + system, + ) + + @staticmethod + def add_extra_box(system): + """ + Add an extra box potential as an obstacle for actin. + """ + if not parameters[f"add_extra_box"]: + return + center = np.array( + [ + parameters[f"extra_box_center_x"], + parameters[f"extra_box_center_y"], + parameters[f"extra_box_center_z"], + ] + ) + size = np.array( + [ + parameters[f"extra_box_size_x"], + parameters[f"extra_box_size_y"], + parameters[f"extra_box_size_z"], + ] + ) + ActinUtil.add_box_potential( + ActinUtil.get_all_actin_particle_types(), + center - 0.5 * size, + size, + ActinUtil.DEFAULT_FORCE_CONSTANT, + system, + ) + @staticmethod def add_dimerize_reaction(system): """ @@ -3246,34 +3405,44 @@ def add_barbed_growth_reaction(system): attach a monomer to the barbed (+) end of a filament. """ for i in ActinUtil.polymer_number_range(): - system.topologies.add_spatial_reaction( - f"Barbed_Growth_ATP1{i}: Actin-Polymer(actin#barbed_{i}) + " - "Actin-Monomer-ATP(actin#free_ATP) -> " - f"Actin-Polymer#GrowingBarbed(actin#{i}--actin#new_ATP)", - rate=parameters["barbed_growth_ATP_rate"], - radius=2 * parameters["actin_radius"] + parameters["reaction_distance"], - ) - system.topologies.add_spatial_reaction( - f"Barbed_Growth_ATP2{i}: Actin-Polymer(actin#barbed_ATP_{i}) + " - "Actin-Monomer-ATP(actin#free_ATP) -> " - f"Actin-Polymer#GrowingBarbed(actin#ATP_{i}--actin#new_ATP)", - rate=parameters["barbed_growth_ATP_rate"], - radius=2 * parameters["actin_radius"] + parameters["reaction_distance"], - ) - system.topologies.add_spatial_reaction( - f"Barbed_Growth_ADP1{i}: Actin-Polymer(actin#barbed_{i}) + " - "Actin-Monomer(actin#free) -> " - f"Actin-Polymer#GrowingBarbed(actin#{i}--actin#new)", - rate=parameters["barbed_growth_ADP_rate"], - radius=2 * parameters["actin_radius"] + parameters["reaction_distance"], - ) - system.topologies.add_spatial_reaction( - f"Barbed_Growth_ADP2{i}: Actin-Polymer(actin#barbed_ATP_{i}) + " - "Actin-Monomer(actin#free) -> " - f"Actin-Polymer#GrowingBarbed(actin#ATP_{i}--actin#new)", - rate=parameters["barbed_growth_ADP_rate"], - radius=2 * parameters["actin_radius"] + parameters["reaction_distance"], - ) + if parameters["barbed_binding_site"]: + bs_number = str(ReaddyUtil.calculate_polymer_number(i, 1, ActinUtil.n_polymer_numbers())) + system.topologies.add_spatial_reaction( + f"Barbed_Growth_BS{i}: Actin-Polymer(binding_site#{i}) + " + "Actin-Monomer-ATP(actin#free_ATP) -> " + f"Actin-Polymer#GrowingBarbed(actin#barbed_ATP_{i}--binding_site#{bs_number})", + rate=parameters["barbed_growth_ATP_rate"], + radius=parameters["binding_site_reaction_distance"], + ) + else: + system.topologies.add_spatial_reaction( + f"Barbed_Growth_ATP2{i}: Actin-Polymer(actin#barbed_ATP_{i}) + " + "Actin-Monomer-ATP(actin#free_ATP) -> " + f"Actin-Polymer#GrowingBarbed(actin#ATP_{i}--actin#new_ATP)", + rate=parameters["barbed_growth_ATP_rate"], + radius=2 * parameters["actin_radius"] + parameters["reaction_distance"], + ) + system.topologies.add_spatial_reaction( + f"Barbed_Growth_ATP1{i}: Actin-Polymer(actin#barbed_{i}) + " + "Actin-Monomer-ATP(actin#free_ATP) -> " + f"Actin-Polymer#GrowingBarbed(actin#{i}--actin#new_ATP)", + rate=parameters["barbed_growth_ATP_rate"], + radius=2 * parameters["actin_radius"] + parameters["reaction_distance"], + ) + system.topologies.add_spatial_reaction( + f"Barbed_Growth_ADP1{i}: Actin-Polymer(actin#barbed_{i}) + " + "Actin-Monomer(actin#free) -> " + f"Actin-Polymer#GrowingBarbed(actin#{i}--actin#new)", + rate=parameters["barbed_growth_ADP_rate"], + radius=2 * parameters["actin_radius"] + parameters["reaction_distance"], + ) + system.topologies.add_spatial_reaction( + f"Barbed_Growth_ADP2{i}: Actin-Polymer(actin#barbed_ATP_{i}) + " + "Actin-Monomer(actin#free) -> " + f"Actin-Polymer#GrowingBarbed(actin#ATP_{i}--actin#new)", + rate=parameters["barbed_growth_ADP_rate"], + radius=2 * parameters["actin_radius"] + parameters["reaction_distance"], + ) system.topologies.add_spatial_reaction( "Branch_Barbed_Growth_ATP1: Actin-Polymer(actin#branch_barbed_1) + " "Actin-Monomer-ATP(actin#free_ATP) -> " @@ -3537,6 +3706,18 @@ def add_translate_reaction(system): rate_function=ReaddyUtil.rate_function_infinity, ) + @staticmethod + def add_position_obstacle_reaction(system): + """ + set the position of the first obstacle particle each timestep. + """ + system.topologies.add_structural_reaction( + "Translate_Obstacle", + topology_type="Obstacle", + reaction_function=ActinUtil.reaction_function_position_obstacle, + rate_function=ReaddyUtil.rate_function_infinity, + ) + @staticmethod def get_position_for_tangent_translation( time_index, displace_stride, monomer_id, monomer_pos, displacement_parameters diff --git a/simularium_readdy_models/common/__init__.py b/simularium_readdy_models/common/__init__.py index 2aa8cfa..b22381d 100644 --- a/simularium_readdy_models/common/__init__.py +++ b/simularium_readdy_models/common/__init__.py @@ -3,3 +3,9 @@ from .particle_data import ParticleData # noqa: F401 from .readdy_util import ReaddyUtil # noqa: F401 from .repeated_timer import RepeatedTimer # noqa: F401 +from .membrane_util import ( + add_membrane_particle_types, + add_membrane_constraints, + init_membrane, + all_membrane_particle_types +) # noqa: F401 diff --git a/simularium_readdy_models/common/membrane_util.py b/simularium_readdy_models/common/membrane_util.py new file mode 100644 index 0000000..0c60e1d --- /dev/null +++ b/simularium_readdy_models/common/membrane_util.py @@ -0,0 +1,282 @@ +#!/usr/bin/env python + +from sys import float_info +import math + +import numpy as np + +from ..common import ReaddyUtil + + +def _leaflet_types(side): + """ + Get all the types for one leaflet side. side = "inner" or "outer" + """ + result = [f"membrane#{side}"] + for n in range(1, 5): + result.append(f"membrane#{side}_edge_{n}") + result.append(f"membrane#{side}_edge_4_1") + result.append(f"membrane#{side}_edge_2_3") + return result + + +def all_membrane_particle_types(): + """ + Get all the membrane particle types. + """ + return _leaflet_types("outer") + _leaflet_types("inner") + + +def add_membrane_particle_types(system, particle_radius, temperature_K, viscosity): + """ + Add particle and topology types for membrane particles + to the ReaDDy system. + """ + diffCoeff = ReaddyUtil.calculate_diffusionCoefficient( + particle_radius, viscosity, temperature_K + ) # nm^2/s + system.topologies.add_type("Membrane") + particle_types = all_membrane_particle_types() + for t in particle_types: + system.add_topology_species(t, diffCoeff) + + +def _add_weak_interaction(types1, types2, force_const, bond_length, depth, cutoff, system): + """ + Adds a weak interaction piecewise harmonic bond to the system + from each type in types1 + to each type in types2 + with force constant force_const + and length bond_length [nm] + and with depth and cutoff. + """ + for t1 in types1: + for t2 in types2: + system.potentials.add_weak_interaction_piecewise_harmonic( + t1, t2, force_const, bond_length, depth, cutoff + ) + + +def _add_box_potential(particle_types, origin, extent, force_constant, system): + """ + Add a box potential to keep the given particle types + inside a box centered at origin with extent. + """ + for particle_type in particle_types: + system.potentials.add_box( + particle_type=particle_type, + force_constant=force_constant, + origin=origin, + extent=extent, + ) + + +def _calculate_lattice(size, particle_radius): + """ + Calculate the x and y lattice coordinates of the membrane + and find the plane dimension with zero size. + """ + plane_dim = -1 + width = [0, 0] + ix = 0 + for dim in range(3): + if size[dim] < float_info.epsilon: + plane_dim = dim + continue + width[ix] = round(size[dim]) + ix += 1 + if plane_dim < 0 or width[1] < float_info.epsilon: + raise Exception("The membrane size must be zero in one and only one dimension.") + result = [ + np.arange(0, width[0], 2. * particle_radius), + np.arange(0, width[1], 2. * particle_radius * np.sqrt(3)) + ] + return result, plane_dim + + +def _calculate_box_potentials(center, size, particle_radius, box_size): + """ + Get the origin and extent for each of the four box potentials + constraining the edges of the membrane. + """ + coords, plane_dim = _calculate_lattice(size, particle_radius) + box_origins = np.zeros((4, 3)) + box_extents = np.zeros((4, 3)) + lattice_dim = 0 + for dim in range(3): + if dim == plane_dim: + box_origins[:, dim] = -0.5 * box_size[dim] + 3 + box_extents[:, dim] = box_size[dim] - 6 + continue + values = coords[lattice_dim] - 0.5 * size[dim] + center[dim] + offset = particle_radius * (1 if lattice_dim < 1 else np.sqrt(3)) + if lattice_dim < 1: + box_origins[0][dim] = values[0] + box_extents[0][dim] = values[-1] - values[0] + box_origins[1][dim] = values[-1] + offset - particle_radius + box_extents[1][dim] = 2 * particle_radius + box_origins[2][dim] = values[0] + offset + box_extents[2][dim] = values[-1] - values[0] + box_origins[3][dim] = values[0] - particle_radius + box_extents[3][dim] = 2 * particle_radius + else: + box_origins[0][dim] = values[0] - particle_radius + box_extents[0][dim] = 2 * particle_radius + box_origins[1][dim] = values[0] + offset + box_extents[1][dim] = values[-1] - values[0] + box_origins[2][dim] = values[-1] + offset - particle_radius + box_extents[2][dim] = 2 * particle_radius + box_origins[3][dim] = values[0] + box_extents[3][dim] = values[-1] - values[0] + lattice_dim += 1 + return box_origins, box_extents + + +def add_membrane_constraints(system, center, size, particle_radius, box_size): + """ + Add bond, angle, and box constraints for membrane particles to the ReaDDy system. + """ + util = ReaddyUtil() + inner_types = _leaflet_types("inner") + outer_types = _leaflet_types("outer") + # weak interaction between particles in the same leaflet + _add_weak_interaction( + inner_types, + inner_types, + force_const=250., bond_length=2. * particle_radius, + depth=7., cutoff=2.5 * 2. * particle_radius, system=system + ) + _add_weak_interaction( + outer_types, + outer_types, + force_const=250., bond_length=2. * particle_radius, + depth=7., cutoff=2.5 * 2. * particle_radius, system=system + ) + # (very weak) bond to pass ReaDDy requirement in order to define edges + util.add_bond( + inner_types, + inner_types, + force_const=1e-10, bond_length=2. * particle_radius, system=system + ) + util.add_bond( + outer_types, + outer_types, + force_const=1e-10, bond_length=2. * particle_radius, system=system + ) + # bonds between pairs of inner and outer particles + util.add_bond( + inner_types, + outer_types, + force_const=250., bond_length=2. * particle_radius, system=system + ) + # angles between inner-outer pairs and their neighbors on the sheet + util.add_angle( + inner_types, + outer_types, + outer_types, + force_const=1000., angle=0.5 * np.pi, system=system + ) + util.add_angle( + inner_types, + inner_types, + outer_types, + force_const=1000., angle=0.5 * np.pi, system=system + ) + # box potentials for edges + box_origins, box_extents = _calculate_box_potentials(center, size, particle_radius, box_size) + corner_suffixes = ["4_1", "2_3"] + for n in range(1, 5): + box_types = [f"membrane#outer_edge_{n}", f"membrane#inner_edge_{n}"] + for suffix in corner_suffixes: + if f"{n}" in suffix: + box_types += [f"membrane#outer_edge_{suffix}", f"membrane#inner_edge_{suffix}"] + break + _add_box_potential( + box_types, + origin=box_origins[n - 1], + extent=box_extents[n - 1], + force_constant=250., system=system + ) + + +def init_membrane(simulation, center, size, particle_radius, box_size): + """ + Add initial membrane particles to the ReaDDy simulation. + """ + coords, plane_dim = _calculate_lattice(size, particle_radius) + cols = coords[0].shape[0] + rows = coords[1].shape[0] + n_lattice_points = cols * (2 * rows) + positions = np.zeros((2 * n_lattice_points, 3)) + types = np.array(2 * n_lattice_points * ["membrane#side-_init_0_0"]) + lattice_dim = 0 + for dim in range(3): + + if dim == plane_dim: + positions[:n_lattice_points, dim] = center[dim] + particle_radius + positions[n_lattice_points:, dim] = center[dim] - particle_radius + continue + + values = coords[lattice_dim] - 0.5 * size[dim] + center[dim] + offset = particle_radius * (1 if lattice_dim < 1 else np.sqrt(3)) + + # positions and types + for i in range(rows): + + p = values if lattice_dim < 1 else values[i] + + start_ix = 2 * i * cols + positions[start_ix:start_ix + cols, dim] = p + if i < 1: + types[start_ix] = "membrane#outer_edge_4_1" + types[start_ix + 1:start_ix + cols] = "membrane#outer_edge_1" + else: + types[start_ix] = "membrane#outer_edge_4" + types[start_ix + 1:start_ix + cols] = "membrane#outer" + + start_ix += n_lattice_points + positions[start_ix:start_ix + cols, dim] = p + if i < 1: + types[start_ix] = "membrane#inner_edge_4_1" + types[start_ix + 1:start_ix + cols] = "membrane#inner_edge_1" + else: + types[start_ix] = "membrane#inner_edge_4" + types[start_ix + 1:start_ix + cols] = "membrane#inner" + + start_ix = (2 * i + 1) * cols + positions[start_ix:start_ix + cols, dim] = p + offset + if i < rows - 1: + types[start_ix:start_ix + cols - 1] = "membrane#outer" + types[start_ix + cols - 1] = "membrane#outer_edge_2" + else: + types[start_ix:start_ix + cols - 1] = "membrane#outer_edge_3" + types[start_ix + cols - 1] = "membrane#outer_edge_2_3" + + start_ix += n_lattice_points + positions[start_ix:start_ix + cols, dim] = p + offset + if i < rows - 1: + types[start_ix:start_ix + cols - 1] = "membrane#inner" + types[start_ix + cols - 1] = "membrane#inner_edge_2" + else: + types[start_ix:start_ix + cols - 1] = "membrane#inner_edge_3" + types[start_ix + cols - 1] = "membrane#inner_edge_2_3" + + lattice_dim += 1 + + membrane = simulation.add_topology("Membrane", types.tolist(), positions) + + # edges + for n in range(n_lattice_points): + other_n = n + n_lattice_points + membrane.get_graph().add_edge(n, other_n) + if n % cols != cols - 1: + membrane.get_graph().add_edge(n, n + 1) + membrane.get_graph().add_edge(other_n, other_n + 1) + if math.ceil((n + 1) / cols) >= 2 * rows: + continue + if (n % (2 * cols) != (2 * cols) - 1): + membrane.get_graph().add_edge(n, n + cols) + membrane.get_graph().add_edge(other_n, other_n + cols) + if (n % (2 * cols) != 0): + membrane.get_graph().add_edge(n, n + cols - 1) + membrane.get_graph().add_edge(other_n, other_n + cols - 1)