Different approaches to computing minimum control energy#
As discussed in Theory, there is more than one way to estimate the minimum control energy associated with a set of state transitions. Here, we compare our standard approach for calculating minimum control energy to one that computes it directly from the controllability Gramian, for many transitions at once. “Why would I want to calculate control energy this way?” I hear you ask. The answer is simple: speed.
The connectome is synthetic, the same kind as in the The protocol workflow tutorial, here with 200 nodes in 20 contiguous systems of 10 nodes each.
import time
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import distance
from nctpy.energies import get_control_inputs, integrate_u, minimum_energy_fast
from nctpy.plotting import reg_plot, set_plotting_params
from nctpy.utils import convert_states_str2int, expand_states, matrix_normalization
set_plotting_params()
def make_connectome(n_nodes=100, n_systems=7, length=0.5, seed=0):
rng = np.random.default_rng(seed)
# node coordinates on a unit sphere, ordered by longitude so that systems are contiguous
xyz = rng.normal(size=(n_nodes, 3))
xyz /= np.linalg.norm(xyz, axis=1, keepdims=True)
xyz = xyz[np.argsort(np.arctan2(xyz[:, 1], xyz[:, 0]))]
names = [chr(ord("A") + i) for i in range(n_systems)]
labels = [names[i] for i, nodes in enumerate(np.array_split(np.arange(n_nodes), n_systems)) for _ in nodes]
# connection probability and weight fall off with distance; within-system pairs are twice as likely
dist = distance.squareform(distance.pdist(xyz, "euclidean"))
decay = np.exp(-dist / length)
p = np.where(np.equal.outer(labels, labels), np.minimum(1, 2 * decay), decay)
edges = np.triu(rng.random((n_nodes, n_nodes)) < p, k=1)
weights = rng.lognormal(0, 1, size=(n_nodes, n_nodes)) * decay
A = np.where(edges, weights, 0)
return A + A.T, dist, labels
A, _, system_labels = make_connectome(n_nodes=200, n_systems=20)
print(A.shape)
(200, 200)
Brain states#
Before we calculate minimum control energy, we need some brain states. We use binary brain states, where each state
comprises a set of regions that are “on” (activity = 1) while the rest of the brain is “off” (activity = 0): one
state per system. states designates which of the 20 states each region belongs to.
states, state_labels = convert_states_str2int(system_labels)
n_nodes = A.shape[0]
n_states = len(state_labels)
print(states)
[ 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 2 2 2 2
2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4
4 4 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 7 7
7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9
9 9 9 9 10 10 10 10 10 10 10 10 10 10 11 11 11 11 11 11 11 11 11 11
12 12 12 12 12 12 12 12 12 12 13 13 13 13 13 13 13 13 13 13 14 14 14 14
14 14 14 14 14 14 15 15 15 15 15 15 15 15 15 15 16 16 16 16 16 16 16 16
16 16 17 17 17 17 17 17 17 17 17 17 18 18 18 18 18 18 18 18 18 18 19 19
19 19 19 19 19 19 19 19]
Minimum energy#
Using these states, we first compute the minimum control energy of every pair of states with our standard
function, get_control_inputs(). Minimum control means the state trajectory is unconstrained
(S all zeros), with every node a control node.
# settings
T = 1 # time horizon
B = np.eye(n_nodes) # all nodes are control nodes
system = "continuous"
A_norm = matrix_normalization(A, system=system)
S = np.zeros((n_nodes, n_nodes)) # x is not constrained
start_time = time.time()
e = np.zeros((n_states, n_states))
for i in np.arange(n_states):
x0 = states == i # ith initial state
for j in np.arange(n_states):
xf = states == j # jth target state
x, u, n_err = get_control_inputs(A_norm=A_norm, T=T, B=B, x0=x0, xf=xf, system=system, S=S)
e[i, j] = np.sum(integrate_u(u))
# integrate_u sums over time points without the time step (dt = 0.001), so divide by 1000 to get the integral
e = e / 1000
print("time elapsed in seconds: {:.2f}".format(time.time() - start_time))
time elapsed in seconds: 15.89
Minimum energy fast#
Now the alternative, minimum_energy_fast(). It needs every transition at once, as a pair of
matrices whose paired columns hold each transition’s initial and target state. expand_states()
builds them from states, for every pair of states:
x0_mat, xf_mat = expand_states(states)
print(x0_mat.shape, xf_mat.shape)
(200, 400) (200, 400)
The rows of x0_mat and xf_mat correspond to the nodes of our system and the columns to the 400
(20 × 20) state transitions. Critically, the two are paired: the same column of both, x0_mat[:, 0] and
xf_mat[:, 0], is the initial and target state of one transition. expand_states only works for binary brain
states; for other states, build x0_mat and xf_mat yourself. Let’s compute energy again:
start_time = time.time()
e_fast = minimum_energy_fast(A_norm=A_norm, T=T, B=B, x0=x0_mat, xf=xf_mat)
e_fast = e_fast.transpose().reshape(n_states, n_states, n_nodes)
e_fast = np.sum(e_fast, axis=2) # sum over nodes
print("time elapsed in seconds: {:.2f}".format(time.time() - start_time))
time elapsed in seconds: 0.52
The fast approach takes a fraction of the time. Is it equivalent?
print(e.shape, e_fast.shape)
print("largest relative difference: {:.1e}".format(np.max(np.abs(e - e_fast) / e_fast)))
f, ax = plt.subplots(1, 1, figsize=(3, 3))
reg_plot(x=e.flatten(), y=e_fast.flatten(), xlabel="minimum energy", ylabel="minimum energy (fast)",
ax=ax, annotate="both", kde=False, regr_line=False)
plt.show()
(20, 20) (20, 20)
largest relative difference: 2.5e-13
The two agree to floating-point precision.
There are caveats. First, minimum_energy_fast() only works for continuous-time systems,
and only for minimum control (no constraint on the state trajectory). Second, it does not output the control
inputs (u), state trajectory (x) or numerical errors (n_err) that get_control_inputs()
provides: you only get node-level energy. Finally, as mentioned above, expand_states() only
works for binary brain states.