MEDYAN.jl
  • Home
  • Tutorials
  • How-To
  • Explanation
  • Reference
  • Versions
    • dev
    • stable

On this page

  • MEDYAN.def_link_type!
    • Tuple{MEDYAN.SysDef}
  • Edit this page
  • Report an issue

MEDYAN.def_link_type!

Tuple{MEDYAN.SysDef}

def_link_type!(s::SysDef; name, description, places, bonds, reactions, param, state)::SysDef

Add a new link type to the system definition. Links represent connections between simulation elements (filament monomers, filament tips, anchors, membrane vertices) that can have associated mechanical bonds and chemical reactions.

Keyword Arguments

  • name::Symbol: Unique identifier for this link type.

  • description::String="": Human-readable description of the link’s purpose.

  • places::Vector{<:Place}: The places this link connects. Each place can be:

    • FilaMonoIdx(): A filament monomer
    • FilaTipIdx(): A filament tip (plus or minus end)
    • Anchor(): A free-floating anchor point (e.g., membrane-bound protein)
    • BallIdx(): A ball (spherical object)
    • MembVertIdx(): A membrane vertex

  • bonds::Vector: Mechanical bonds between places. Each bond is a NamedTuple with:

    • bond: Bond type (e.g., DistanceRestraint(), custom Bond subtype)
    • input::Tuple{Int...}: Which places the bond connects
    • param::NamedTuple: Fixed parameters (e.g., spring constant k)
    • state::NamedTuple: Mutable state (e.g., rest length L0)
    • enabled::Bool=true: Whether bond is active by default

  • reactions::Vector{Vector}: Reactions for each place. reactions[i] is a Vector of reactions for places[i]. Trailing places with no reactions can be omitted (i.e. length(reactions) ≤ length(places)). Each reaction is a NamedTuple with:

    • Required fields:

      • name::Symbol: Unique identifier for this reaction within the place.
      • affect!: Callback that executes the reaction. Signature: (c::Context; link::Link, chem_voxel::Int, reaction_id, place_idx, kwargs...) -> Int. link is a handle identifying the specific link instance — pass it to get_state, update_link!, link2tags, etc. Return a status code for debugging purposes. For fila_cutoff reactions, also receives place::FilaMonoIdx (a randomly chosen nearby filament monomer within the cutoff distance). For ball_cutoff reactions, also receives place::BallIdx (a randomly chosen nearby ball within the cutoff distance).

    • Optional fields:

      • rate = Returns(1.0): Function (c::Context; link::Link, link_data, place_idx, link_state::NamedTuple, kwargs...) -> Float64 returning the state-dependent rate factor. link_state is the current mutable state of this link instance (as defined by the state keyword). For multi-site binding, return the number of available sites (e.g., MAX_ARP23 - link_state.num_arp23). link_data is an internal object — pass it through to link2tags or get_link_mechanics for a small performance gain over the 1-argument forms. See “When rates are reevaluated” below for important constraints on what this function should depend on.
      • base_rate::Float64 = 1.0: Rate constant multiplier. Units depend on invvolumepower: 1/s for unimolecular (invvolumepower=0), nm³/s for bimolecular (invvolumepower=1). For example, a diffusion-limited on-rate k_D = 4πDR ≈ 1.57×10⁹ nm³/s for D=25 μm²/s, R=5 nm. Separating large rate constants from small integer multipliers in rate avoids overflow.
      • invvolumepower::Int64 = 0: Volume scaling. Set to 1 for bimolecular (nm³/s), 0 for unimolecular (1/s).
      • reactants_extra::String = "": Additional reactant species whose counts multiply the propensity (e.g., "diffusing.ARP23"). Uses the same "+"-separated format as def_reaction! (e.g., "diffusing.A + diffusing.B").
      • fila_cutoff::Tuple{Symbol, Float64} = nothing: If set, e.g., (:actin, 50.0), the propensity is multiplied by the number of nearby filament monomers of that type within the cutoff distance (nm). When the reaction fires, affect! receives one such monomer as place::FilaMonoIdx. Some monomers are invisible to this search until after a mechanics cycle — see “Nearby count caching and mechanics” below.
      • fila_tip_cutoff::Tuple{Symbol, Bool, Float64} = nothing: If set, e.g., (:actin, true, 50.0), the propensity is multiplied by the number of nearby filament tips of that type within the cutoff distance (nm). The Bool is is_minus_end: true for minus-end tips, false for plus-end tips. When the reaction fires, affect! receives one such tip as place::FilaTipIdx. Mutually exclusive with fila_cutoff and ball_cutoff. Some tips are invisible to this search until after a mechanics cycle — see “Nearby count caching and mechanics” below.
      • ball_cutoff::Float64 = nothing: If set, e.g., 50.0, the propensity is multiplied by the number of nearby balls whose centers are within the cutoff distance (nm). When the reaction fires, affect! receives one such ball as place::BallIdx. Only balls that have been through a mechanics cycle (is_minimized) are visible to this search. Mutually exclusive with fila_cutoff and fila_tip_cutoff.
      • enabled::Bool = true: Whether reaction is active by default.

  • param::NamedTuple=(;): Fixed parameters for the link type.

  • state::NamedTuple=(;): Initial mutable state (e.g., (num_bound=Int32(0),)).

When rates are reevaluated

Although the rate callback receives the full Context and could in principle read any simulation state, the cached rate value is only incrementally recalculated for a specific link when one of the following events occurs during a chemistry step:

  1. update_link! is called on this link — any change to state, places, bond_states, bond_enabled, or reaction_enabled triggers a rate recalculation for all reactions on this link.

  2. The chemical state (get_chem_state) of an attached place changes — specifically:

    • A filament monomer state changes (e.g., via polymerization, depolymerization, severing, or update_fila_mono_state!): rates are recalculated for links attached to the changed monomer and its ±1 neighbors (and attached tips).
    • An anchor is moved via update_anchor!.
    • A ball is modified via update_ball!.

The entire cache is fully rebuilt at the start of each chemistry phase after mechanics — so all rates are recomputed from scratch using the current Context state after mechanics.

This means: within a chemistry step, if the rate function reads something from the Context other than link_state or the chem state of its own places (see get_chem_state), that dependency will not be tracked and the cached rate will become stale until the next full cache rebuild. For example:

  • ✅ Safe: Reading link_state fields (e.g., link_state.num_arp23). Updated immediately when update_link! changes state.
  • ✅ Safe: Reading get_chem_state of attached places, or calling get_link_mechanics / is_minimized on the same link (these change only after mechanics, which triggers a full cache rebuild before the next chemistry step).
  • ⚠️ Stale within a chemistry step: Reading diffusing species counts, other links’ states, or the chem state of places that this link is not attached to. These values may change during chemistry without triggering a rate update for this link.

Workaround — rejection sampling: If a reaction rate genuinely depends on state that rate cannot track (e.g., the state of a nearby monomer not attached to this link), have rate return an upper bound of the true rate. Then in affect!, read the current state and accept the reaction with probability true_rate / upper_bound (rejecting by returning 0 otherwise). This preserves correct stochastic kinetics despite the stale cache.

Nearby count caching and mechanics (fila_cutoff / fila_tip_cutoff / ball_cutoff reactions)

When a reaction has fila_cutoff set, its total propensity is base_rate × rate(...) × nearby_monomer_count × [reactants_extra counts] / V^invvolumepower. When a reaction has fila_tip_cutoff set, the propensity is the same but with nearby_tip_count in place of nearby_monomer_count. When a reaction has ball_cutoff set, the propensity is the same but with nearby_ball_count in place of nearby_monomer_count. These counts are cached alongside the rate factor and have important update semantics:

Snapshotted positions

1. After each mechanics cycle, monomer and tip positions are snapshotted.
   These snapshotted positions do not change during chemistry.
2. The nearby count is computed using the link place's **current** position
   (i.e. the `Anchor`, `FilaMonoIdx`, etc. whose reaction has the cutoff)
   against the snapshotted positions. The count is recomputed whenever the
   link place moves or a new link is created.
3. Monomers that have not yet been through a mechanics cycle (from
   polymerization, filament creation, or severing) are not in the snapshot.
4. Similarly, a tip is visible only if it existed at the last mechanics
   snapshot. Polymerization and depolymerization move an existing tip
   but do not create a new one, so the tip stays visible at its
   snapshotted position.
5. A ball is visible to `ball_cutoff` only when
   [`is_minimized`](@ref) is `true`. Changing `position` or `radius`
   via [`update_ball!`](@ref) resets `is_minimized` to `false`.

Example: Membrane-bound NPF cluster

const NPF_MAX_ARP23 = Int32(3)  # maximum ARP2/3 that can bind per NPF

MEDYAN.def_link_type!(s;
    name=:npf_anchor,
    description="Nucleation promoting factor anchored to membrane",
    places=[Anchor()],
    state=(;
        num_arp23 = Int32(0),  # bound ARP2/3 count
        num_actin = Int32(0),  # bound actin count
    ),
    reactions=[
        [  # reactions for places[1] (the Anchor)
            # Bimolecular binding: NPF + diffusing ARP2/3 → NPF·ARP2/3
            (;
                name = :arp23_bind,
                affect! = (c; link, chem_voxel, kwargs...) -> let
                    link_state = get_state(c, link)
                    update_link!(c, link; state=(num_arp23=link_state.num_arp23 + Int32(1),))
                    add_diffusing_count!(c; species=:ARP23, chem_voxel, amount=-1)
                    1  # success
                end,
                rate = (c; link_state, kwargs...) -> NPF_MAX_ARP23 - link_state.num_arp23,
                base_rate = 5E7,  # nm³/s (bimolecular on-rate)
                reactants_extra = "diffusing.ARP23",
                invvolumepower = 1,
            ),
            # Unimolecular unbinding: NPF·ARP2/3 → NPF + ARP2/3
            (;
                name = :arp23_unbind,
                affect! = (c; link, chem_voxel, kwargs...) -> let
                    link_state = get_state(c, link)
                    update_link!(c, link; state=(num_arp23=link_state.num_arp23 - Int32(1),))
                    add_diffusing_count!(c; species=:ARP23, chem_voxel, amount=+1)
                    1  # success
                end,
                rate = (c; link_state, kwargs...) -> link_state.num_arp23,
                base_rate = 0.01,  # 1/s per bound ARP2/3
            ),
        ],
    ],
)
  • Edit this page
  • Report an issue