Solving ODEEdges with algebraic variables without mass matrix

[Antomek]

Hello all,

First I would like to thank you for the work being done on NetworkDynamics.jl. It looks very promising and potentially very useful!

I had a question on how to model heterogeneous networks. I hope this is the right place to ask it.
Specifically, suppose I have two types of vertices: type A and type B.
A vertex can be connected to either a vertex of the same type of the other type.

The network is described by a directed graph, and the connection is specified by the vertex from which it originates.
So essentially we have two types of edges: ones originating from type A vertices, and those originating from type B vertices.

Ideally I would like to do something of the lines of (in the ODE describing a vertex):

for e in e_d
    if source_vertex == typeA
        dv[1] += connection_A(v, e)
    elseif source_vertex == typeB
        dv[2] += connection_B(v, e)
    end
end

but at the moment I have no idea how to achieve something like this.

Any help is appreciated!

[lindnemi]

Hello Antomek,

you could define two types of edges and pass them as an edge list to the network_dynamics constructor. For that you need to get an ordered list of the edges in the graph via LightGraphs.jl collect(edges(graph)). Then you create a list of edge functions that has an edgeA! function if the corresponding edge orginates at a vertex A, resp edgeB for a vertex B. You construct the network via network_dynamics(vertex_list, edge_list, graph). The coupling term (summation) in the vertex can then simply iterate over every element in e_d. Does that work for your problem?

Since we are always excited when our package find a new use case, I would be very curious to learn in which domain your system arises. Could you point me to any references?

[Antomek]

Hello Lindnemi,

Thank you for the reply.
I am playing around with using your package to simulate neurodynamics.
One principle in neurophysiology, Dale's principle, states that neurons release the same set of neurotransmitters at their synapses.

If we model the synapses as the edges, this means that a neuron (or vertex) has the same electric effect on the downstream vertices with which it is connected.
The desire to have different edges stems from the existence of different neuron types with different neurotransmitters, and by the above principle, the dynamic effect of an edge is determined by the vertex from which it originates.

The upshot is that ideally I would like to mix an ODEEdge and StaticEdge or something along those lines.
More concretely: suppose that I would like to do in the ODEEdge:

function edgeA!(de, e, v_s, v_d, p, t)
    de[1] = f(v_s, e[1]) # ODE describing synaptic conductance depending on presynaptic voltage
    e[2] = g(e[1], v_d) # Electric current which gets added to downstream neuron's voltage derivative
end

and then be able to do something in the vertex along the lines of

for e in e_d
    dv[1] += e[2]
end

I'm sure a workaround is possible but I cannot see it for the moment.
Hopefully I've somewhat clarified what I'm trying to achieve here.
In spirit it really is quite similar to the FitzHughNagumo model that you have in the examples!

[Lima-Lima]

Hi Antomek,

Great application! I work with @lindnemi on the package. We've discussed this exact sort of situation before and are aware that it is useful for models in many different domains. You're correct to observe that at the moment there isn't a clean way to make it work with the package. It's in our planned features.

In fact, we do have a prototype in place that implements what you describe. The example similar to yours is here from a working branch on my fork. Ignore the comment about DAE, that results from an earlier discussion. Take note that this branch isn't stable or complete, but the dynamics you describe are possible. We are still discussing how to finish up the interface to this feature.

Are you interested in taking the proof-of-concept to a polished pull request for this feature, or should we let you know when we have completed it? It will move up the feature queue for me now that you've opened this issue.

@lindnemi , I think we should rename or merge this issue to represent the "plasticity" issue, because the heterogeneous modeling in the original question you answered nicely.

[Antomek]

Hey Lima-Lima,

Thank you, that example does look promising!
I am still a GitHub novice, so will take a look into making a pull request myself.
First I'll take a look at the prototype.
Thanks for the help!

[lindnemi]

There is already another way implemented to make the above version of edgeA work. You have to use a static edge with a mass matrix. Internally this will solve the equation 0 = M[2,2] * de[2] = e[2] - g(e[1], v_d) and thus satisfiy the constraint described in the edge function.

function edgeA!(de, e, v_s, v_d, p, t)
    de[1] = f(v_s, e[1]) # ODE describing synaptic conductance depending on presynaptic voltage
    e[2] = g(e[1], v_d) # Electric current which gets added to downstream neuron's voltage derivative
end

M = zeros(2,2)
M[1,1] = 1

nd_edgeA! = ODEEdge(f! = edgeA!, dim= 2, mass_matrix = M)

However this has two drawbacks:

1. You have to use a solver which is capable of handling mass matrices (e.g. Rodas4, Rosenbrock23) and those are usually slower than explicit algorithms. [They may be sped up by letting ModelingToolkit symbolically derive a Jacobian matrix for the problem]a

2. You have to choose the initial conditions in a way that satisfies the algebraic constraints. Otherwise solve will have to find a valid inital conditions before starting the integration and this might take ages or fail altogether.

With these issues in mind these types of mixed edges work fine, especially for smaller systems where optimized performance is not an issue. However the prototype by @Lima-Lima will eventually sidestep these issues altogether and any help getting this to a pull request is greatly appreciated :)

[Antomek]

Is there any update on doing this?
I keep running into this problem of having heterogeneous edges, and needing the class or "type" of the connection in the vertex dynamics.

This suggestion:

you could define two types of edges and pass them as an edge list to the network_dynamics constructor. For that you need to get an ordered list of the edges in the graph via LightGraphs.jl collect(edges(graph)). Then you create a list of edge functions that has an edgeA! function if the corresponding edge orginates at a vertex A, resp edgeB for a vertex B.

unfortunately doesn't work if you, e.g. have different values stored in the edges.
For example, if you're looking to do something along the lines of:

for e in edges
    if class(e) == :classA
        dv[1] += e[1]
    elseif class(e) == :classB
        dv[2] += e[2]
    end
end

Are there ways of solving this that I am overlooking?

[lindnemi]

Unfortunately no update on this. Next year we will have more funding to advance NetworkDynamics, at the moment the package is mostly in maintenance mode. If you want to open a PR on this problem, that would be very welcome. Anyways I thought of the following workaround:

Basically you want to try to unify the addition, i.e., use

for e in edges
    dv[1] += e[1]
    dv[2] += e[2]
end

There are several ways to achieve that. If you do not need to keep the edge variables, you can overwrite those which are not used in the addition, e.g.,

function classB(e, v_s, v_d, K, t)
    e[1] = f(...)
    e[2] = g(...)
    e[1] = 0.
end

If you do indeed need them, you can increase the edge dimension

function classB(e, v_s, v_d, K, t)
    e[1] = f(...)
    e[2] = g(...)
    e[3] = 0.
    e[4] = e[2]
end

in combination with

for e in edges
    dv[1] += e[3]
    dv[2] += e[4]
end

Does this help with your problem?