Threshold Linear Networks

These examples are based on the paper Graph Rules for Recurrent Neural Network Dynamics by Carina Curto and Katherine Morrison.

using LinearAlgebra
using SparseArrays
using Catlab
using Catlab.Graphs
using OrdinaryDiffEq
using NonlinearSolve
using AlgebraicDynamics
using AlgebraicDynamics.ThresholdLinear
using Test
using Plots
using Catlab.Graphics

Triangle Graph

trig = @acset Graph begin
  V = 3
  E = 3
  src = [1,2,3]
  tgt = [2,3,1]
end

Graph {V:3, E:3}

E	src	tgt
1	1	2
2	2	3
3	3	1

tri = TLNetwork(CTLNetwork(trig))
@test tri.W' == [ 0.0   -0.75  -1.5;
          -1.5    0.0   -0.75;
          -0.75  -1.5    0.0;
      ]
prob = ODEProblem(tri, [1,1/2,1/4], [0.0,40])
soln = solve(prob, alg=Tsit5())
plot(soln)

Test case that u1 = 0 still converges to attractor with support 123

prob = ODEProblem(tri, [0,1/2,1/4], [0.0,40])
soln = solve(prob, alg=Tsit5())
plot(soln)

support(soln)

3-element Vector{Int64}:
 1
 2
 3

Symmetric Edge Graph

We want to build a test case for non-full support. We look to the symmetric edge graph.

barbell = CTLNetwork(@acset Graph begin
  V = 2
  E = 2
  src = [1,2]
  tgt = [2,1]
end)

CTLNetwork{Float64}(Graph:
  V = 1:2
  E = 1:2
  src : E → V = [1, 2]
  tgt : E → V = [2, 1], LegalParameters{Float64}(0.25, 0.5, 1.0))

From the theory, we know that adding an isolated vertex shouldn't affect the fixed points. The isolated vertex has no activators, so it will decay to zero, and it activates nothing, so it can't change those fixed points.

bb3 = CTLNetwork(apex(coproduct([barbell.G, Graph(1)])))
prob = ODEProblem(bb3, [1/2, 1/4, 1/8], (0.0,100))
soln = solve(prob, Tsit5())
@test support(soln, 1e-14) == [1,2]
plot(soln)

Embedded subgraphs should have their fixed points persist. We add the extra vertex adjacent to both vertices in the barbell.

barbell_plus = CTLNetwork(@acset Graph begin
  V = 3
  E = 4
  src = [1,2,3,3]
  tgt = [2,1,1,2]
end)
TLNetwork(barbell_plus).W

3×3 Matrix{Float64}:
  0.0   -0.75  -0.75
 -0.75   0.0   -0.75
 -1.5   -1.5    0.0

Notice how symmetric this W matrix is.

prob = ODEProblem(barbell_plus, [1/2, 1/4, 1/8], (0.0,100))
soln = solve(prob, Tsit5())
@test support(soln, 1e-14) == [1,2]
plot(soln)

The categorical coproduct of graphs gives you the direct sum of their adjacency matrices. When you direct sum graphs, you should cartesian product their fixed points. Two disjoint symmetric edges have fixed points on each.

barbell_pair = CTLNetwork(apex(coproduct(barbell.G, barbell.G)))
draw(barbell_pair.G)

We first look for the attractor supported on the 3rd and 4th vertex (σ=[3,4]).

prob = ODEProblem(barbell_pair, [1/2, 1/4, 2/3, 4/5], (0,50))
soln = solve(prob, Tsit5())
@test support(soln, 1e-12)== [3,4]
plot(soln)

With a different initial condition, we find another attractor this one corresponding to σ=[1,2].

prob = ODEProblem(barbell_pair, [1/2, 1/4, 2/50, 1/40], (0,50))
soln = solve(prob, Tsit5())
@test support(soln, 1e-12)== [1,2]
plot(soln)

The graphs that have interesting structure have a lot of symmetry so we try to make one with a product. Categorical products have symmetry because they work like the cartesian product of sets.

bt = apex(product(add_reflexives(barbell.G), trig))
tln = CTLNetwork(bt)
draw(tln.G)

Now we look for an attractor.

prob = ODEProblem(tln, [1/2, 1/4, 2/3, 4/5, 1/10, 5/6], (0,150))
soln = solve(prob, Tsit5())
@test support(soln, 1e-5) == [1,2,4,5,6]
plot(soln)

We can find another attractor by zeroing out some variables in the initial condition.

prob = ODEProblem(tln, [0, 1/4, 2/3, 0, 1/10, 0], (0,150))
soln = solve(prob, Tsit5())
plot(soln)

Notice that even when you have only a singleton in the initial condition, you don't get a singleton support in the attractor.

prob = ODEProblem(tln, [0, 0, 2/3, 0, 0, 0], (0,150))
soln = solve(prob, Tsit5())
# @test support(soln, 1e-5)  == [2,3,6]
plot(soln)

Because of symmetry in the model, we can pick out a different attractor.

prob = ODEProblem(tln, [0, 2/3, 0, 0, 0, 0], (0,150))
soln = solve(prob, Tsit5())
# @test support(soln, 1e-5) == [2,3,4]
plot(soln)

Using Nonlinear Solvers to find fixed points

NonlinearSolvers.jl lets us define the steady state of our system as our fixed point. We want unstable fixed points, so we can't use the DynamicSS problem type provided by SciML. We have to use traditional root finders rather than an evolve to equilibrium approach.

prob = NonlinearProblem(tln, [0, 2/3, 0, 0, 0, 0])
fp = solve(prob)
fp.u

6-element Vector{Float64}:
 0.14285714285714282
 0.14285714285714293
 0.1428571428571429
 0.1428571428571428
 0.1428571428571429
 0.1428571428571428

Once we compute the fixed point, we can plug it in to the dynamics and simulate. This finds the corresponding oscillatory attractor due to the numerical perturbations. We could converge to this attractor faster by adding a perturbation to fp.u.

prob = ODEProblem(tln, fp.u, (0,150))
soln = solve(prob, Tsit5())
plot(soln)

Because root finders are only guaranteed to find local minima of the residual, we start at a different initial guess and find a different attractor.

prob = NonlinearProblem(tln, [1/3, 0, 1/2, 1/2, 1/2, 1/2])
fp = solve(prob)

retcode: Success
u: 6-element Vector{Float64}:
 0.30769230769230776
 1.198286922568818e-18
 0.17582417582415233
 0.17582417582419932
 0.10989010989007968
 0.10989010989014007

We can plug in the fixed point and find the oscillatory attractor.

prob = ODEProblem(tln, fp.u, (0,150))
soln = solve(prob, Tsit5())
plot(soln)

Induced Subgraphs Preserve Attractors

When you take an induced subgraph. You can restrict the dynamics onto that subgraph.

g = induced_subgraph(bt, [1,2,3])
tln = TLNetwork(CTLNetwork(g))
prob = NonlinearProblem(tln, [1/3, 0, 1/2])
fp = solve(prob)
support(fp.u)
prob = ODEProblem(tln, fp.u, (0,150))
soln = solve(prob, Tsit5())
plot(soln)

The following two functions automate the analysis that we did above

Restrict to a subgraph
Solve for a fixed point in the subgraph
Plug that solution in to the dynamics of the full system
Solve those dynamics and plot

AlgebraicDynamics.ThresholdLinear.restriction_fixed_point — Function

restriction_fixed_point(G::AbstractGraph, V::AbstractVector{Int}, parameters=DEFAULT_PARAMETERS)

Restrict a graph to the induced subgraph given by V and then solve for its fixed point support.

Returns the fixed point of G corresponding to nonzeros in V padded with zeros for the vertices that aren't in V.

source

function restriction_simulation(G, V, tspan=(0,150.0), parameters=DEFAULT_PARAMETERS)
  tln = CTLNetwork(G, parameters)
  σ, u₀ = restriction_fixed_point(G, V)
  @show u₀
  @show σ
  prob = ODEProblem(tln, u₀, tspan)
  soln = solve(prob, Tsit5())
  plt = plot(soln)
  return σ, soln, plt
end

restriction_simulation (generic function with 3 methods)

Mining patterns in our product graph

Let's take a look at our graph again

draw(bt)

We can try finding an attractor from the triangle 2,4,6

σ₀, soln, plt = restriction_simulation(bt, [2,4,6]);
plt

Here we can look at a tiny subgraph.

σ₀, soln, plt = restriction_simulation(bt, [1,2]);
plt

There is a bigger support on 2,4,5,6

σ₀, soln, plt = restriction_simulation(bt, [2,4,5,6]);
plt

Trying 1,2,4

σ₀, soln, plt = restriction_simulation(bt, [1, 2,4]);
plt

Trying 1,2,3,5

σ₀, soln, plt = restriction_simulation(bt, [1,2,3,5]);
plt

Library Reference

You can access these functions from the module AlgebraicDynamics.ThresholdLinear.

AlgebraicDynamics.ThresholdLinear.DEFAULT_PARAMETERS — Constant

DEFAULT_PARAMETERS = (ϵ=0.25, δ=0.5, θ=1.0)

source

AlgebraicDynamics.ThresholdLinear.CTLNetwork — Type

CTLNetwork{F}

G::Graph
parameters::LegalParameters{F}

Stores a combinatorial linear threshold network as a graph and the parameters. The single argument constructor uses the DEFAULT_PARAMETERS.

source

AlgebraicDynamics.ThresholdLinear.LegalParameters — Type

LegalParameters{F}

The parameters of a CTLN model are:

epsilon
delta
theta

The constructor enforces the constraint that ϵ < δ/(δ+1).

source

AlgebraicDynamics.ThresholdLinear.TLNetwork — Type

 TLNetwork{T}

W::AbstractMatrix{T}
b::AbstractVector{T}

Stores a Threshold Linear Network as a Matrix of weights and Vector of biases. You can construct these from a Graph with CTLNetwork.

source

AlgebraicDynamics.ThresholdLinear.TLNetwork — Method

TLNetwork(g::CTLNetwork)

Convert a CTLN to a TLN by realizing the weights an biases for that Graph.

source

AlgebraicDynamics.ThresholdLinear.dynamics — Method

dynamics(tln::TLNetwork)

Construct the vector field for the TLN. You can pass this to an ODEProblem, or as a Continuous Resource Sharer.

source

AlgebraicDynamics.ThresholdLinear.indicator — Method

indicator(g::AbstractGraph, σ::Vector{Int})

Get the indicator function of a subset with respect to a graph. This should probably be a FinFunction.

source

AlgebraicDynamics.ThresholdLinear.nldynamics — Method

nldynamics(tln::TLNetwork)

Construct the root finding problem for the TLN. You can pass this to an NonlinearProblem to find the steady states of the network.

source

AlgebraicDynamics.ThresholdLinear.restriction_fixed_point — Function

restriction_fixed_point(G::AbstractGraph, V::AbstractVector{Int}, parameters=DEFAULT_PARAMETERS)

Restrict a graph to the induced subgraph given by V and then solve for its fixed point support.

Returns the fixed point of G corresponding to nonzeros in V padded with zeros for the vertices that aren't in V.

source