wsphillips / conductor.jl Goto Github PK

Network is no longer working after the updates to ConductanceSystem and CompartmentSystem. It should be fixed by writing a NetworkSystem that similar implements the ModelingToolkit system interface.

The default constructor for NetworkSystem should support taking a Metagraph or similar datastrcuture capable of describing network topology, neuron types, synapse types and weights.

For now we will create a directly coupled system of equations. Future features may include spike-triggered discrete synaptic events.

Instead of verbosely writing out an electrode stimulus (e.g. holding current, pulse trains), we should be able to offer standard constructors that let users specify stimuli with the same style as pulse train and waveform generators used in electrophysiology. We can of course support other types of stimuli, but I'd like to have at least a core standard offering that works like a Master-8/9

There should be an easier and more readable way to specify gating variables, channels, compartments, etc.

Ref #8 (comment) for beginnings of this.

We should support 1-dimensional branched cable models. This boils down to creating a conversion from a CompartmentSystem to PDESystem in ModelingToolkit.jl. We should be able to convert any CompartmentSystem with a cylindrical Geometry to a discretized PDE, given some delta X and boundary conditions.

MulticompartmentSystem is a thin layer over CompartmentSystem that allows connecting compartments together. Combined with the above, it should be enough to get going with branched cable models.

To support these systems we also need:

• Spatially dependent parameters (e.g. gbar for dendritic channel types may have a gradient as a function of x)
• Augment synapses and stimuli so that we can target specific locations (e.g. given a cylindrical span of length 50µm, place a synapse/electrode at x = 40µm etc). This could be done by inserting synapse/stimuli terms to specific compartments post discretization or by using a boolean/dirac delta type term to turn them on at a set location.
• Support for reading in external data (namely SWC files; to be tracked in a separate issue. I have an MWE of this already)

(Some) possible approaches to implementation:

1. Discretize compartments ourselves. Divide a CompartmentSystem into many individual compartments, apply modifications (e.g. synapses, parameter values, etc) and then broadcast convert.(ODESystem, comps)
2. Convert from CompartmentSystem to PDESystem and use MethodOfLines.jl to discretize. This is likely more robust but I'm not sure how to setup branched compartments with this approach.

Finally, it will be interesting to see how performant out-of-the-box code generation is compared to existing domain-specific optimizations (e.g. Hines matrix + many examples of problem-specific decomposition strategies for parallelization etc)

Take advantage of #57 to make problem setup easier for downstream applications

There's an @info macro in the Simulation function which gives details about the model and shows up after solving the simulation. For example:
Model CaS with 3 equations
States (4):
g(t)
h(t) [defaults to (1.0 + exp(9.67742 + 0.16129Vₘ(t)))^-1]
m(t) [defaults to (1.0 + exp(-4.07407 - (0.123457Vₘ(t))))^-1]
Vₘ(t)
Parameters (1):
gbar [defaults to 4.0]
...

However, this message keeps showing up in the Julia REPL after running any command (even after writing something completely unrelated to the neuron problem, like x=1). Any idea for why is this happening?

There are scaling issues with larger simulations (~1000+ neurons) during model construction and code generation. Models are slow to compile or may crash the compiler.

An example of the above can be seen in demo/pinksy_network.jl. If n_neurons is incremented to 500+ neurons, the compiler will hang for a very long time on the subsequent call to solve, possibly crashing before reaching the solution. On my Intel i5-12600k, it took ~40 minutes to solve for 500 neurons.

If the problem is scaled even further, to ~1000 neurons, code generation hard fails at the problem construction step due to inference failure related to numerous callbacks. Type inference hits a stack overflow because of the extremely high number of discrete callbacks (one for each neuron). Unfortunately, this is unavoidable due to the need to scalarize neurons combined with a lack of vectorized callback constructor for discrete events, analogous to VectorContinuousCallback

From what I can tell these are two separate problems but both arise from scalarizing components (neurons) for ModelingToolkit.jl

Until the (hopefully imminent) release of improved symbolic IR in ModelingToolkit, I recommend avoiding very large simulations and, where possible, using continuously integrated callback events for synapse models in larger systems. Continuous callbacks are more computationally expensive during model execution, but should scale better in terms of compilation since we can use VectorContinuousCallback instead of a very long CallbackSet

Unit checking in ModelingToolkit is now available and checks occur automatically for symbols with unit metadata. These new features should replace or augment most of the existing units handling currently in Conductor.jl

Is direct support for ustrip all that's needed for

Just as with channels, compartments should also be implemented as a custom system type.

CompartmentSystem needs to implement conversion to both ODESystem and PDESystem, but we can start with just ODESystem.

Quick notes:

• Instead of putting equations, states, params etc in single containers, we should divide up storage across multiple fields and lazily evaluate methods like equations(::CompartmentSystem. This will let us organize and selectively modify certain things more easily. For instance, it simplifies the logic behind pushing and popping additional membrane currents if we have a field like currents::Vector that we use to generate the current sum term for D(Vm) on demand.
• We should have a better building block, analogous to GatingVariable, that can be used for defining state transformations (i.e. a replacement for AuxConversion)
• In the ideal case, components in a compartment should validate whether their inputs are satisfied by the parent compartment. This should include querying whether an input type could be imputed from other available states/transformations in cases where it might not be explicitly provided by another component. MTK has the notion of inputs/outputs started, but not fully implemented afaict. So we might want to refrain from getting too carried away on this point.
• More notes to be added as edits to this post.

Instead of starting with an adjacency matrix, build MulticompartmentSystem on a tree data structure (potentially AbstractTrees.jl interface). It will be more robust for equation building, and it should eliminate the need to depend on parent references directly in CompartmentSystem.

Later, during lowering, we will likely transform to an adjacency matrix for PDE Cable equation modeling. Again, having a tree to start from will make this easier since a Hines matrix is a depth first search on a dendritic tree.

Use the the Pinsky & Rinzel model as a template to write a test set for MultiCompartmentSystem

Test MultiCompartmentSystem with more compartments--make a demo and test set for the complete Traub 1991 19-compartment model and Davison 2000 4-compartment model.

Hi @wsphillips !

I tried to run the STGneuron file in the demo folder, with a freshly downloaded Conductor, and get this error:

WARNING: could not import Conductor.Cationic into Main
ERROR: LoadError: MethodError: no method matching MembranePotential(::Quantity{Int64, 𝐋² 𝐌 𝐈⁻¹ 𝐓⁻³, Unitful.FreeUnits{(mV,), 𝐋² 𝐌 𝐈⁻¹ 𝐓⁻³, nothing}})
Closest candidates are:
  MembranePotential() at ~/.julia/packages/Conductor/idY1O/src/Conductor.jl:95

Stacktrace:
 [1] top-level scope
   @ ~/.julia/dev/Conductor/demo/prinz_kinetics.jl:6
in expression starting at /Users/andrearamirez/.julia/dev/Conductor/demo/prinz_kinetics.jl:6

These are the packages in my environment

Status `~/.julia/dev/Conductor/demo/Project.toml`
  [d194a933] Conductor v0.0.2
  [615f187c] IfElse v0.1.1
  [961ee093] ModelingToolkit v6.4.9
  [1dea7af3] OrdinaryDiffEq v6.7.1
  [91a5bcdd] Plots v1.27.1
  [295af30f] Revise v3.3.3
  [1986cc42] Unitful v1.11.0

and I'm using Julia Version 1.7.0.

The same error shows up when I try to run STGnetwork.jl and hodgkinhuxley.jl. These demos should work straight from the src code, right?

Thanks for your help!

Hi @wsphillips,

I'm trying to model a neuron with conductances regulated by an integral control rule that uses the calcium concentration, [Ca]. See https://www.cell.com/neuron/pdf/S0896-6273(14)00292-X.pdf
The new equations that I want to add are eqs (3) in the paper:
D(g_i) = 1/tau_g * (mRNA_i - g_i)
D(mRNA_i) = 1/tau_i * (Ca_{target} - [Ca])

I think that ConductanceSystem could contain the ODESystem with the equations for D(g) and D(mRNA) as well as D(gating variables) but I don't know how to implement these changes. What would be the best way to incorporate them into the code ?

Thanks!

Certain type parameters are recruiting the type system for "switch"-like usage rather than optimized code paths. An example is ion types, which can just be something simple like an enum

We should be able to connect together assemblies of compartments. These types of neuron models follow the form of those similar to:

Pinsk & Rinzel 1994 -- two compartment model

Davison et al 2000 -- four compartments

and range up in size to (for instance) the classic 19-compartment Traub et al 1992 model.

These types of models could be implemented by defining a function that specifies equations for the conductance between individual compartments and returns a new system using compose from MTK. Additional tooling to handle compartment-specific synaptic conductances and stimuli will also be necessary (e.g. post-synaptic connection targeted to Neuron Y, compartment Y)

We should allow construction and manipulation of neuronal networks as a graph object instead of explicit adjacency lists (e.g. Vector{Synapse}). Ideally we should support multilayer graphs that are compatible with the Graphs.jl interface

Would be great if we had functions to easily connect neurons directly, by specifying the dynamics of the synapse and (in the case of phenomenological models) specifying to which compartment the synapse goes.

For Networks this could be extended so that one can use an adjacency matrix / probability of connection presence, or loop over all the neurons and specify the connection between individual pairs.

From parent issue #66

Add a cached value for the the time stamp of the last spike event for each neuron so that we can support time-based refractory periods in addition to the current simple check for positive-going threshold crossings. A spike would then occur only when all of the following conditions are met:

1. V > V_threshold
2. V_previous < V_threshold
3. t - t_lastspike > refractory_period

When a spike occurs then the associated effect should additionally reset the value of t_lastspike

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Thank you for this great package and bringing ModelingToolkit to the world of neuronal modeling! I was interested in simple threshold based neuron models (QIF, Izhikevich), so I started out with the LIF demo in order to adopt it later on. But it seems sth changed in the meantime, causing the following error with the current version @main and Julia 1.8:

ERROR: MethodError: no method matching CompartmentSystem(::Conductor.LIF; name=:stim_neuron)

And in general, how would you describe the current state of this package? E.g. would you already welcome contributions for other models to build up a zoo or is there a lot of change in the overall structure to be expected, so this would only cause overhead in adjustments all the time?

Allow "threshold only" gating of synaptic events (and therefore postsynaptic currents). Currently, spike events are limited to positive-going threshold-crossing events. This behavior enforces a refractory period until the voltage of the presynaptic root compartment (i.e. the soma) goes back below the threshold voltage.

The kinetics of some models depend on a continuous/repetitive suprathreshold activation. To accomodate these, we need an option to repeatedly trigger suprathreshold events.

Splitting out from #16

Allow specifying distributions for variables in lieu of or in combination with default values.

Fallback behavior should be to randomly draw from a distribution if running a simulation with no default values specified. Eventually exploit this information to modify lowering step when targeting Gen.jl

Currently spike events are hard coded to occur when:

1. Somatic voltage potential exceeds 10mV.
2. If the compartment is refractory, a spike occurs only if the membrane potential was below threshold on the previous time step. For continuous spike events, this is implicit since the affect is definied for positive-going threshold events.

Both of the above are placeholders. We should ultimately support changing the spike threshold and/or defining arbitrary conditions to be met for a spike event. We should also support setting refractory periods based on elapsed time or other state-based conditions.

Changing the raw threshold for spikes is easy--we can extend the existing spike detection methods with an extra optional argument. At minimum, a time-based refractory period requires caching the time stamp of the last spike. You would then define an event as (pseudocode) V > Vthold && Vprev < Vthold && t - t_lastspike > refractory

Support for arbitrary state conditions could be achieved if a user defines their own functors for spike detection. At minimum, we can enable this by defining abstract supertypes for users to overload. However, this approach is burdensome for the user since it requires working on a lower level, using the vanilla DiffEq callbacks API, instead of the ModelingToolkit.jl symbolic callback interface. I will need to think about this last component some more and likely address it as a separate issue.

From parent issue #66

Spike thresholds are hard coded to occur at 10mV membrane potential. Add an extra argument during construction of synaptic models or compartment models that lets users adjust the value of the threshold.

We should start using the new compose and extend provided by ModelingToolkit.jl in lieu of the older systems=[sys1,sys2] syntax. Using extend we can more simply express things like adding applied/electrode current to existing neuron models, as well as plugging in synaptic currents.

We should be able to easily take the dynamics of a single neuron, and then construct a whole group of these neurons.

We need to figure out how to address the individual neurons in each group (e.g. through indexing), and introduce some variations within the group when desired.

This issue is a key part of the work necessary for solving #61

Provide a keyword argument for defining additional affects that modify the state of presynaptic CompartmentSystem following a spike-triggered event. The default use case is for voltage resetting (e.g. in LIF neurons). This will let us move away from using "baked-in" model-specific constructors.

ModelingToolkit.jl offers some additional metadata tags that we can use to more easily setup optimization. We want to be able to:

• Mark certain symbols as constants
• Mark symbols as "tunable" for inclusion in the state vector for optimization packages.
• Bounds and distributions for symbols
• Mark symbols as inputs/outputs

In models where spikes are determined by a simple threshold crossing event, we should optionally support a refractory period of arbitrary length. It should also be a hook for other state manipulations. Practically, this could be handled by chaining an additional callback affect to spike events.

Using the version of Conductor.jl installed from the REPL (0.0.4, Julia 1.7.3), the voltage step response generated by the Simulation function is a factor 1000x faster than expected after generating the neuron model via the CompartmentSystem constructor. This situation can also be generated in the most recent hodgkinhuxley.jl demo script by removing the active channels. Minimal example notebook showing this below:

https://gist.github.com/rpgowers/b1fff319f8832736287e8788aedb9e11

Take advantage of the callback features in ModelingToolkit.jl to create synaptic currents that reset based on presynaptic threshold crossing.