On the control of recurrent neural networks using constant inputs

Cyprien Tamekue; Ruiqi Chen; ShiNung Ching

doi:10.1109/tac.2025.3615934

. Author manuscript; available in PMC: 2026 Mar 24.

Published in final edited form as: IEEE Trans Automat Contr. 2025 Sep 30;71(3):1737–1752. doi: 10.1109/tac.2025.3615934

On the control of recurrent neural networks using constant inputs

Cyprien Tamekue ¹, Ruiqi Chen ², ShiNung Ching ³

PMCID: PMC13008370 NIHMSID: NIHMS2152725 PMID: 41878724

Abstract

This paper investigates the controllability of a broad class of recurrent neural networks widely used in theoretical neuroscience, including models of large-scale human brain dynamics. Motivated by emerging applications in non-invasive neurostimulation such as transcranial direct current stimulation (tDCS), we study the control synthesis of these networks using constant and piecewise constant inputs. The neural model considered is a continuous-time Hopfield-type system with nonlinear activation functions and arbitrary input matrices representing inter-regional brain interactions. Our main contribution is the formulation and solution of a control synthesis problem for such nonlinear systems using specific solution representations. These representations yield explicit algebraic conditions for synthesizing constant and piecewise constant controls that solve a two-point boundary value problem in state space up to higher-order corrections with respect to the time horizon. In particular, the input is constructed to satisfy a tractable small-time algebraic relation involving the Jacobian of the nonlinear drift, ensuring that the synthesis reduces to verifying conditions on the system matrices. For canonical input matrices that directly actuate $k$ nodes, this implies that the reachable set (with constant inputs) of a given initial state is an affine subspace whose dimension equals the input rank and whose basis can be computed efficiently using a thin QR factorization. Numerical simulations illustrate the theoretical results and demonstrate the effectiveness of the proposed synthesis in guiding the design of brain stimulation protocols for therapeutic and cognitive applications.

Index Terms—: Control Synthesis, Nonlinear Control, Recurrent Neural Networks, Neurostimulation, Transcranial Direct Current Stimulation

I. Introduction

Non-invasive neurostimulation techniques, such as transcranial magnetic stimulation (TMS) and transcranial electrical stimulation (tES), are increasingly being used to modulate brain activity in order to achieve desired cognitive outcomes [4], [15], [22]. Conceptualizing the brain as a controlled dynamical system [24], [34] offers a powerful framework for identifying key brain regions and networks [26] involved in specific cognitive functions and understanding how they can be modulated using targeted stimulation.

In this regard, there is strong potential at the intersection of network control theory and clinical and basic neuroscience. For instance, there have been successes in the application of linear network control theory to optimize TMS, administered with a constant input over short time intervals [25]. Such applications target the control of specific brain regions based on structural network connectivity derived from diffusion spectrum imaging (DTI). Similarly, at an analysis level, linear network control theory has been employed in human neuroimaging studies, again by leveraging DTI-parameterized models [16]. In this latter study, the authors postulated how specific brain regions may be suitable targets for control based on advantageous controllability properties relative to their potential actuation. A complementary perspective views brain stimulation—particularly deep brain stimulation—as a form of vibrational control that stabilizes desired neural dynamics through high-frequency inputs, a notion that has been formalized and validated in recent theoretical work [27], [28].

While these earlier studies [16], [25] have provided valuable insights by using linear network control theory and structural connectivity data to understand brain dynamics, they are limited in key respects. Foremost, by virtue of linearity, they are only assured to provide local characterizations. In [25], it is shown that these local analyses can offer some predictive power over how stimulation induces functional changes in nonlinear models. However, the changes considered are in terms of low-dimensional correlation metrics between brain areas (network nodes), as opposed to specific configurations in state space. In contrast, our goal is to enable formal control synthesis to induce arbitrary network states in the presence of full nonlinearity, which is likely essential for understanding how different brain dynamics mediate cognitive function.

Indeed, the need for more biologically realistic, nonlinear models of large-scale brain dynamics has been appreciated [5] as a precursor for the application of control theory to cognitive neuroscience. Among the paradigms in this regard are whole-brain Mesoscale Individualized NeuroDynamic (MINDy) models [8], [32]. Our previous works have demonstrated the validity and utility of MINDy models as a generative tool for understanding the relationship between individualized neural architectures and neural dynamics in both resting-state [32], [33] and cognitive task contexts [31]. The MINDy models derived from single-subject resting-state brain imaging data operate at a macroscale level, where each node represents a distinct brain region. A given model captures the temporal evolution of brain activity and the non-linear interactions between hundreds of brain regions, taking the form of a nonlinear dynamical system of continuous-time Hopfield-type recurrent neural networks [18]. All (intrinsic) parameters–the decay and connectivity matrices and the transfer function parameters–are directly estimated from brain activity time series so that the resulting models predict future brain activity, providing a more accurate and biologically plausible representation of large-scale brain dynamics.

A MINDy model belongs to the more general class of Hopfield-type recurrent neural networks, which will constitute the focus of our study herein. Specifically, we consider a model of interconnected neural masses associated with $d$ regions, described by a set of nonlinear differential equations that represent the dynamics of each region’s state over time. The dynamic for the $i$ -th region in this network evolves as

\frac{d x_{i}}{d t} (t) = - α_{i} x_{i} (t) + \sum_{j = 1}^{d} W_{i j} f_{j} (x_{j} (t)) + \sum_{j = 1}^{k} b_{i j} u_{j} x_{i} (0) = x_{i}^{0}

(1)

where $k \in N^{*}$ satisfies $k \leq d, x_{i} (t)$ denotes the $i$ -th region’s state at time $t \geq 0, α_{i} > 0$ is a positive parameter reflecting the rate at which the current state of the $i$ -th region decays, $u_{j}$ are constant control inputs, $b_{i j}$ represent coupling coefficients between control inputs and neurons, $W_{i j}$ is a real constant that weights the connection from the $j$ -th region to the $i$ -th region, and $f_{j}$ is the activation function of the $j$ -th region.

The study of neural networks, such as (1), has long been a topic of great interest for researchers seeking to understand complex dynamical behaviors and their associated control mechanisms, [3], [13], [18], [36], [40], [41]. Among the various properties investigated, controllability–the ability to steer a network from one state to another within a finite time–stands out for its significant theoretical and practical implications. As highlighted earlier, in fields such as neuroengineering and brain stimulation, particularly transcranial electrical stimulation (tES) with direct current (stimulation via a constant input), the operative issue is to determine a constant input that can drive the state of the network (1) to a desired target state–within the state space–over a short time horizon (Figure 1). Such capability would present a promising avenue for developing stimulation protocols tailored to achieve specific therapeutic outcomes [12], such as altering the excitability of the motor cortex [4].

Fig. 1. — We consider the control of network models of the general form (1). We are motivated by emerging applications in neurostimulation involving the manipulation of brain networks by means of exogenous input.

In previous works such as [13], [36], [40], [41], controllability or stability properties of equations as (1) have been considered in the full-actuated configuration when $d = k$ and $B = (b_{i j}) = Id$ or in the specific cases of $α_{i} = 0$ for all $i$ , and $B$ belonging to an open dense subset of $ℳ_{d, k} (R)$ .

In the current study, we do not impose a priori such a restriction on the input matrix $B \in ℳ_{d, k} (R)$ . The assumptions on the activation functions are relaxed to those commonly used in the literature in studying neural networks. The general consideration on $B$ and the complexity arising from the inherent nonlinear dynamics present significant challenges, and systematic solutions for the control synthesis problem are not readily available in the current literature. To address this challenge, we derive an explicit algebraic representation for the hidden state $x (t)$ and the associated constant input that solves the two-point boundary value problem in state space up to higher-order terms. This yields an explicit small-time algebraic condition for the constant input, showing that the required control can be determined by inverting an expression involving the drift’s Jacobian and the desired change in state, up to higher-order corrections with respect to the time horizon. From this derivation, we then provide a step function (piecewise constant control) for large periods of time. Although tDCS typically relies on constant or slowly modulated inputs, step constant controls provide a theoretically grounded and numerically efficient extension of constant-input synthesis over longer time horizons. Note that, in the case of linear activation functions, the synthesized constant control remains valid for both short and long time horizons.

For input matrices that directly actuate $k$ nodes, we show that the reachable set with constant inputs is an affine subspace whose dimension equals the rank of $B$ , and its basis can be computed explicitly via a thin QR factorization of the associated constraint matrix. This characterization applies to both linear and nonlinear activation functions, clarifying how the local structure of the input matrix shapes the short-horizon reachable set.

The remainder of the paper is organized as follows. In the next section, we introduce the general notations used throughout the paper. Section II presents our modeling framework, reformulating equation (1) into a form more suitable for controllability analysis and establishing key results concerning solution representations. Section III addresses the control synthesis problem and is divided into several parts: Section III-A focuses on the linear activation case; Section III-B develops control synthesis strategies for general nonlinear activations that appear to be the proper generalization of the linear case; and Section III-C outlines a taxonomy of synthesis approaches derived from dual representations of the solution, clarifying the conceptual origin of the forward nominal-state and backward nominal-state strategies. Section III-D then provides practical considerations by analyzing the structure of reachable sets associated with the proposed step-function control syntheses. Section IV offers broader discussion and possible extensions. Section V presents numerical illustrations of the theory. Finally, Section VI summarizes the main results and discusses directions for future research. Technical proofs of some results are deferred to the Appendix.

Notation 1. In this paper, $Z$ stands to the set of integers, $N^{*}$ denote the set of positive integers and $N = N^{*} \cup {0}$ . For all $a, b \in N$ , we denote by $⟦a, b⟧ = N \cap [a, b]$ . For all $d, k \in N^{*}$ , we denote by $ℳ_{d, k} (R)$ the set of $d \times k$ matrices with real coefficients. $R^{d}$ denotes the $d$ -dimensional real column vector and $ℒ (R^{d})$ the space of linear maps on $R^{d}$ that we identify, in the usual way, with $ℳ_{d} (R)$ the set of squared matrices of order $d$ with real coefficients. For all $(x, y) \in {(R^{d})}^{2}$ , we denote by $| x |$ the Euclidean norm of $x$ and by $⟨ x, y ⟩$ the scalar product of $x$ and $y$ . We denote the identity matrix in $ℳ_{d} (R)$ by Id, and for every matrix $M \in ℳ_{d, k} (R)$ , $‖ M ‖ = sup \{| M x | |x \in R^{k}, | x | = 1\}$ denotes the spectral norm of $M$ , and $M^{†}$ denotes the Moore-Penrose pseudo-inverse of $M$ . For a symmetric matrix $M \in ℳ_{d} (R)$ , $λ_{min} (M)$ and $λ_{max} (M)$ denote respectively its smaller and largest eigenvalues. For ease of notation, a diagonal matrix $Γ = diag (γ_{1}, γ_{2}, \dots, γ_{d}) \in ℳ_{d} (R)$ , will simply denoted as $Γ = (γ_{i} δ_{i, j})$ , where $δ_{i, j}$ is the usual Kronecker symbol. If $A \in ℳ_{d} (R)$ , then $e^{A} \in ℳ_{d} (R)$ and $σ (A)$ are, respectively, the exponential and the set of eigenvalues of $A$ , viz.

e^{A} = \sum_{n = 0}^{\infty} \frac{A^{n}}{n!}, σ (A) = \{λ \in C∣ ker (A - λ Id) \neq 0\} .

(2)

We recall that $𝒪$ stands to the big $𝒪$ -notation, and $y (t) \underset{t ~ 0}{=} 𝒪 (x (t))$ means $y (t)$ is of the order of $x (t)$ . Finally, we use $\dot{x} (t) = \frac{d x}{d t} (t)$ for the total derivative of $x$ w.r.t. $t$ .

II. Representation of solutions

To study the controllability of (1), it is convenient to recast it in abstract form. Introduce the (nonlinear) map

N (x) = - D x + W f (x) \forall x \in R^{d} .

(3)

Then, (1) recasts as

\dot{x} (t) = N (x (t)) + B u, x (0) = x^{0}

(4)

where $x = {(x_{1}, x_{2}, \dots, x_{d})}^{⊤} \in R^{d}$ is the state’s vector, $x^{0} = {(x_{1}^{0}, \dots, x_{d}^{0})}^{⊤} \in R^{d}$ is the initial state, $u = {(u_{1}, \dots, u_{k})}^{⊤} \in R^{k}$ is a constant control, $B = (b_{i j}) \in ℳ_{d, k} (R)$ is the input matrix, $D = (α_{i} δ_{i, j}) \in ℳ_{d} (R)$ is the decay matrix, $W = (W_{i j}) \in ℳ_{d} (R)$ is the connectivity matrix and $f : R^{d} \to R^{d}$ given by $f (x) = {(f_{1} (x_{1}), f_{2} (x_{2}), \dots, f_{d} (x_{d}))}^{⊤}$ is the firing rate function of the network.

Notation 2. We set $A = - D + W$ , where $D$ is the decay matrix and $W$ is the connectivity matrix in (3).

Assumptions II.1. Unless otherwise stated, we assume that for every $i \in ⟦1, d⟧$ , the activation function $f_{i} : R \to R$ is a $C^{2}$ , globally Lipschitz function on $R$ with a bounded second derivative, and satisfies $f_{i} (0) = 0$ . The latter is without loss of generality since we may always take for every $s \in R, {\tilde{f}}_{i} (s) = f_{i} (s) - f_{i} (0)$ and set $\tilde{B u} = B u + W f (0)$ as the new input to equation (4). Finally, for the sake of simplicity in the presentation, we also assume that ${‖f_{i}^{'}‖}_{\infty} = 1$ . Indeed, as long as ${‖f_{i}^{'}‖}_{\infty} \neq 0$ , we can always define an activation function ${\tilde{f}}_{i} (s) = f_{i} (λ s)$ with $λ = 1 / {‖f_{i}^{'}‖}_{\infty}$ and $s \in R$ .

Throughout the following, we will refer to the following fixed parameters:

Γ = λ_{min} (D) - ‖W‖, Λ = λ_{max} (D) + ‖W‖,

Λ_{1} = ‖ W ‖ {‖f^{″}‖}_{\infty} and {‖f^{″}‖}_{\infty} = max_{i \in 〚 1, d 〛} {‖f_{i}^{″}‖}_{\infty} .

Recall that the vector field $N$ is globally Lipschitz on $R^{d}$ (see, for instance, Lemma A.1). Denoting by $ϕ_{t} : R^{d} \to R^{d}$ its flow at $t \in R$ , it follows that the family $\{ϕ_{t} ∣ t \in R\}$ is a one parametric subgroup of $Diff (R^{d})$ , the group of diffeomorphism of $R^{d}$ . Note that for any $x^{0} \in R^{d}, ϕ_{t} (x^{0})$ is the unique solution of (4) when the control $u = 0$ , namely

\frac{d}{d t} ϕ_{t} (x^{0}) = N (ϕ_{t} (x^{0})), ϕ_{0} (x^{0}) = x^{0} .

(5)

Let $ψ_{t}$ denote the inverse of $ϕ_{t}$ , i.e., $ψ_{t} = {(ϕ_{t})}^{- 1} = ϕ_{- t}$ . Then, for every $x^{0} \in R^{d}$ , it holds

\frac{d}{d t} ψ_{t} (x^{0}) = - N (ψ_{t} (x^{0})), ψ_{0} (x^{0}) = x^{0} .

In particular, the map $(t, x) \in R \times R^{d} \mapsto ϕ_{t} (x) \in R^{d}$ is of class $C^{2}$ (see, for instance, [17, Chapter 15, Theorem 1]). Moreover, for a fixed $t \in R$ , if we use $D ϕ_{t} (x)$ to denote the differential of $ϕ_{t}$ , then $D ϕ_{t} (x)$ is a well-defined invertible matrix, and it holds

{[D ϕ_{t} (x)]}^{- 1} = D ψ_{t} (ϕ_{t} (x)), \forall x \in R^{d} .

(6)

The classical theory of ordinary differential equations (ODE) can be invoked to justify the existence and uniqueness of solutions of (4) (see, for instance, [6, Chapter 2]). In this work, we represent the system’s state trajectory in a form reminiscent of that of linear time-invariant systems as used in [38, Theorem 4.4]. This representation emerged from the chronological calculus framework introduced by Agrachev and Gamkrelidze in the late 1970s for solving nonautonomous ODEs on finite-dimensional manifolds [2]. For a comprehensive treatment of this theory, see also [1, Chapter 6], and [21], [30] for applications in geometric control theory. This representation is noteworthy because, in the specific case of linear activation functions, it naturally aligns with the variation of the constants formula. From this perspective, the proposed representation and approach are a proper generalization of the variation of constants formula to the case of models of the form (4).

The first result of this paper is about solution representation, and it states as follows. The proof is presented in Section B-A.

Theorem II.2. Let $T > 0$ and $(x^{0}, u) \in R^{d} \times R^{k}$ . The solution $x \in C^{3} ([0, T]; R^{d})$ of (4) can be expressed as

x (t) = ϕ_{t} (x^{0} + \int_{0}^{t} D ψ_{s} (x (s)) B u d s), \forall t \in [0, T] .

(7)

When the dynamics are linear, this representation reduces to the familiar form.

Corollary II.3. Let $T > 0$ and $(x^{0}, u) \in R^{d} \times R^{k}$ . If $f_{i} (s) = s$ , the solution $x \in C^{\infty} ([0, T]; R^{d})$ to (4) reads

x (t) = e^{t A} x^{0} + \int_{0}^{t} e^{(t - s) A} B u d s \forall t \in [0, T] .

(8)

Proof. In this case, $N (x) = A x$ for every $x \in R^{d}$ where $A$ is introduced in Notation 2. In this case, $ϕ_{t} (x) = e^{t A} x$ and $ψ_{t} (x) = e^{- t A} x$ for every $t \in R$ . It follows that $D ψ_{t} (x (t)) = e^{- t A}$ and the result then follows by (7). □

We conclude this section by presenting a second representation of the solution to (4). It can be viewed as a natural extension of the linear case given in Corollary II.3, particularly at the final time horizon $T$ . The proof follows similar lines to that of Theorem II.2.

Theorem II.4. Let $T > 0$ and $(x^{0}, u) \in R^{d} \times R^{k}$ . Then the solution $x \in C^{3} ([0, T]; R^{d})$ of (4) can be expressed as

x (t) = ψ_{T - t} (ϕ_{T} (x^{0}) + \int_{0}^{t} D ϕ_{T - s} (x (s)) B u d s)

(9)

for all $t \in [0, T]$ .

Remark II.5. Although Theorems II. 4 and II. 2 are stated for constant inputs $u \in R^{k}$ , they remain valid for time-dependent inputs $u \in L^{1} ([0, T]; R^{k})$ . Note that this is a specific case of [37, Theorem 3.2]. In this case, the corresponding solution $x (\cdot)$ to (4) belongs to $C^{0} ([0, T]; R^{d})$ . This observation is particularly relevant, as we will later consider step constant controls—i.e., piecewise constant controls defined on $[0, T]$ .

III. Controllability of the neural networks with STEP FUNCTION

In this section, we analyze the controllability of (4) under constant inputs over short time horizons and extend the approach to step-function controls for longer-horizon transfers in the fully nonlinear setting. When the activation functions are linear, the synthesized constant control remains valid for both short and long time horizons.

Recall, e.g., from [39, Remark 4.2.3] that a step function on $[0, T]$ is a piecewise constant function defined over $[0, T]$ , i.e, constant on each interval of a partition of $[0, T]$ . Then, a constant function corresponds to a special case of a step function with a single constant value across the entire interval.

Let us introduce the following.

Definition 3 (Step controllability). Let $T > 0$ . System (4) is step controllable over the time interval $[0, T]$ if, for all $x^{0}, x^{1} \in R^{d}$ , there exists a step function $u$ on $[0, T]$ such that the solution of (4) with $x (0) = x^{0}$ satisfies $x (T) = x^{1}$ .

A. The case of linear activation functions

To gain insight into the controllability of the nonlinear system (4) under constant control, it is instructive to first analyze the corresponding linearized model—a strategy that has been clearly articulated and effectively motivated in prior work, notably by [19].

Consider first the linear neural network

\dot{x} (t) = A x (t) + B u, x (0) = x^{0} \in R^{d},

(10)

with $A = - D + W$ (corresponding to $f_{i} (s) = s$ for all $i$ ) and constant input $u \in R^{k}$ . System (10) is controllable on $[0, T]$ if and only if the Kalman rank condition holds. Equivalently, the invertibility of the controllability Gramian also guarantees controllability over $[0, T]$ . Moreover, when these conditions are met, one can synthesize a time-varying control $\underline{u} \in L^{\infty} ((0, T); R^{k})$ of minimal $L^{2}$ -energy that steers (10) from any $x^{0}$ to any $x^{1}$ in time $T$ (see, e.g., [10], [39]).

However, this synthesis does not systematically extend to constant control capable of achieving the same objective. If the system is controllable, and we assume the existence of a constant input $\tilde{u} \in R^{k}$ that steers the system from $x^{0} \in R^{d}$ to $x^{1} \in R^{d}$ over the interval $[0, T]$ , then from (8) and $x (T) = x^{1}$ ,

\int_{0}^{T} e^{(T - s) A} d s B \tilde{u} = x^{1} - e^{T A} x^{0}

so that left multiplying by $A$ and using the identity

A \int_{0}^{T} e^{(T - s) A} d s = \int_{0}^{T} (- \frac{d}{d s} e^{(T - s) A}) d s = e^{T A} - Id

(11)

we find

(e^{T A} - Id) B \tilde{u} = A (x^{1} - e^{T A} x^{0}) .

(12)

However, (12) can be ill-posed, as shown by the following counter-example.

Remark III.1. Consider the case where $k = d = 2, T = 2 π / ω, B = Id$ , and

A = [\begin{matrix} 0 & - ω \\ ω & 0 \end{matrix}], e^{T A} = [\begin{matrix} cos (2 π) & - sin (2 π) \\ sin (2 π) & cos (2 π) \end{matrix}] = Id .

The linear system (10), equipped with time-varying control, is fully actuated and hence controllable, as guaranteed by the Kalman rank condition. However, any constant control $\tilde{u} \in R^{2}$ must satisfy (12). Since $e^{T A} - Id = 0$ , it follows that no constant control $\tilde{u} \in R^{2}$ can steer $x^{0} = 0$ to a nonzero target state $x^{1} \neq 0$ over the interval $[0, T]$ . Notably, $T = 2 π / ω$ can be made arbitrarily small by choosing $ω ≫ 2 π$ .

On the other hand, using the Gramian-based synthesis, one finds that $W_{c}^{- 1} = Id / T$ , and

\underline{u} (t) = \frac{1}{T} e^{(T - t) A^{⊤}} x^{1}, t \in [0, T]

is the time-varying control of minimal $L^{2}$ -norm that steers $x^{0} = 0$ to $x^{1} \neq 0$ . In particular, $‖ \underline{u} ‖_{L^{2}}^{2} = ⟨W_{c}^{- 1} x^{1}, x^{1}⟩ = {|x^{1}|}^{2} / T$ . Taking $x^{1} = (\sqrt{T / ω}, 0)^{⊤}$ yields $‖ \underline{u} ‖_{L^{2}}^{2} = 1 / ω$ , showing that if $ω ≫ 1$ , then $\underline{u}$ achieves the transfer from $x^{0} = 0$ to $x^{1} = (\sqrt{T / ω}, 0)^{⊤}$ in $T > 0$ with small effort.

Remark III.1 suggests that the controllability of a linear system does not imply the existence of a constant control that solves the control objective.

Proposition III.2. Let $T > 0$ and $(x^{0}, x^{1}) \in {(R^{d})}^{2}$ . Assume that for any $ℓ \in Z, i \frac{2 π ℓ}{T} \notin σ (A)$ . Then, the solution $x (\cdot)$ to (10) satisfies $x (T) = x^{1}$ , if and only if, $u \in R^{k}$ solves

B u = {(e^{T A} - Id)}^{- 1} A (x^{1} - e^{T A} x^{0}) .

(13)

The proof of the following result is immediate.

Corollary III.3. Under hypotheses of Proposition III.2, (13) has at least one solution $u \in R^{k}$ , if and only if,

{(e^{T A} - Id)}^{- 1} A (x^{1} - e^{T A} x^{0}) \in Im B .

In this case, $u \in R^{k}$ given by

u = B^{†} {(e^{T A} - Id)}^{- 1} A (x^{1} - e^{T A} x^{0})

is the corresponding least-norm constant control.

Proof of Proposition III.2. Under the assumption on the eigenvalues of $A$ , the matrices $A$ and $e^{T A}$ – Id are invertible. By Corollary II.3, the solution $x \in C^{\infty} ([0, T]; R^{d})$ of (10) with a constant control $u \in R^{k}$ is given by (8). Let us show that $x (T) = x^{1}$ is equivalent to (13). If $x (T) = x^{1}$ , then one obtains immediately that $u \in R^{k}$ solves (13) by invoking the same arguments used to derive (12). Conversely, if (13) is satisfied, one deduces from (11) and (8) that

A x (T) = A e^{T A} x^{0} + (e^{T A} - Id) B u = A x^{1}

so that $x (T) = x^{1}$ since $A$ is invertible. □

Remark III.4. If $‖ W ‖ < λ_{min} (D)$ , then the matrix $e^{T A} - Id$ is invertible for all $T > 0$ , as ensured by Lemma C.2 and the Neumann series expansion. Moreover, for all $x \in R^{d}$ , we have

| A x | \geq | D x | - | W x | \geq (λ_{min} (D) - ‖ W ‖) | x |

which implies that $A$ is invertible as well.

Therefore, in Proposition III.2, the spectral assumptions become relevant only when $‖ W ‖ \geq λ_{min} (D)$ .

Remark III.5. Remark III. 1 and Proposition III. 2 highlight a key distinction between constant and time-varying control synthesis in linear systems.

While time-varying control synthesis depends solely on the structural properties of the pair ( $A, B$ )–as captured by the invertibility of the controllability Gramian—constant control synthesis is more restrictive. In particular, it depends not only on the matrices $A$ and $B$ , but also crucially on the time horizon $T > 0$ , the initial and target states $x^{0}$ and $x^{1}$ , and on the spectrum of $A$ . As shown in Remark III.1, even a fully actuated system may fail to connect certain states under constant input when $e^{T A} = Id$ , rendering (12), then (13) unsolvable.

B. The case of nonlinear activation functions

Proposition III.2 highlights that controllability under constant control depends intricately on the time horizon $T$ , the spectrum of $D N (\cdot)$ , and the initial and target states $x^{0}$ and $x^{1}$ .

Based on the solution representations of (4) given in Theorems II.2 and II.4, we note that the latter naturally extends the linear case, particularly at the final time horizon $T$ . Nonetheless, both formulations provide valid foundations for synthesizing controls that achieve the desired state transfer. The key distinction lies in their temporal structure and directionality, as well as in the specific conditions required to guarantee their applicability.

1) Forward nominal-state synthesis: Building on the solution representation of (4) given in Theorem II.4, we aim to synthesize step-function controls that achieve the control objective. Unlike the linear case, the control is first derived in an implicit form due to the nonlinear nature of the system. Then, under the assumption of a small time horizon, we provide an explicit expression that is an accurate approximation of the step-function controls that solve the control objective.

We begin by stating a key result on which the constant control synthesis is based. The proof is identical to that presented in Section C-A.

Proposition III.6. Let $T > 0$ and $(x_{0}, u) \in R^{d} \times R^{k}$ . The solution $x (\cdot)$ of (4) can be expanded at time $t = T$ as

x (T) = \sum_{n = 0}^{\infty} (- 1)^{n} \frac{T^{n + 1} D N {(ϕ_{T} (x^{0}))}^{n}}{(n + 1)!} D ϕ_{T} (x^{0}) B u + ϕ_{T} (x^{0}) - φ_{u} (T) B u .

(14)

Here $φ_{u} (T) = κ_{u} (T) + η_{u} (T)$ ,

κ_{u} (T) = \sum_{n = 1}^{\infty} \int_{0}^{T} \frac{(t - T)^{n + 1}}{(n + 1)!} [\frac{d}{d t} Z_{T - t}^{n}] Q_{T - t} d t,

(15)

η_{u} (T) = \sum_{n = 1}^{\infty} \int_{0}^{T} \frac{(t - T)^{n} Z_{T - t}^{n - 1}}{n!} D^{2} ϕ_{T - t} (x (t)) \dot{x} (t) d t

(16)

where $Q_{T - t} = D ϕ_{T - t} (x (t))$ , $Z_{T - t} = D N (ϕ_{T - t} (x (t)))$ and $D^{2} ϕ_{T - t} (x (t))$ is the second derivative of $ϕ_{T - t}$ at $x (t)$ .

The implicit constant control synthesis of this section is therefore presented in the following theorem.

Theorem III.7. Let $T > 0$ and $(x^{0}, x^{1}) \in {(R^{d})}^{2}$ . Assume that for any $ℓ \in Z, i \frac{2 π ℓ}{T} \notin σ (D N (ϕ_{T} (x^{0})))$ . Then, the solution $x (\cdot)$ to (4) satisfies $x (T) = x^{1}$ , if and only if, $u \in R^{k}$ solves

[Id - 𝒜_{T} (x^{0}) φ_{u} (T)] B u = 𝒜_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0})) .

(17)

Here $φ_{u} (T) \in ℳ_{d} (R)$ is defined in Proposition III.6, and letting $U_{T} (x^{0}) = D N (ϕ_{T} (x^{0}))$ , one has

𝒜_{T} (x^{0}) = D ψ_{T} (ϕ_{T} (x^{0})) {(Id - e^{- T U_{T} (x^{0})})}^{- 1} U_{T} (x^{0}) .

(18)

Proof. Let us show that $x (T) = x^{1}$ is equivalent to (17). If $x (T) = x^{1}$ , then left multiplying (14) by $U_{T} (x^{0})$ , using (2), and

\sum_{n = 0}^{\infty} \frac{{(- T U_{T} (x^{0}))}^{n + 1}}{(n + 1)!} = \sum_{n = 1}^{\infty} \frac{{(- T U_{T} (x^{0}))}^{n}}{n!} = e^{- T U_{T} (x^{0})} - Id

yields

U_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0})) = - U_{T} (x^{0}) φ_{u} (T) B u + (Id - e^{- T U_{T} (x^{0})}) D ϕ_{T} (x^{0}) B u .

(19)

Since ${[D ϕ_{T} (x^{0})]}^{- 1} = D ψ_{T} (ϕ_{T} (x^{0}))$ , we left multiply (19) respectively by $({Id - e^{- T U_{T} (x^{0})})}^{- 1}$ and $D ψ_{T} (ϕ_{T} (x^{0}))$ to get

𝒜_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0})) = [Id - 𝒜_{T} (x^{0}) ζ_{u} (T)] B u

which is exactly (17) where $𝒜_{T} (x^{0})$ is defined by (18). Conversely, assume that $u \in R^{k}$ solves (17). Then, one finds

𝒜_{T} (x^{0}) φ_{u} (T) B u = B u - 𝒜_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0})) .

(20)

One the other hand, from (14), one deduces that

𝒜_{T} (x^{0}) φ_{u} (T) B u = B u - 𝒜_{T} (x^{0}) (x (T) - ϕ_{T} (x^{0})) .

(21)

Identifying (20) with (21), one gets

𝒜_{T} (x^{0}) (x (T) - x^{1}) = 0 .

(22)

Since $𝒜_{T} (x^{0})$ is invertible, it follows from (22) that $x (T) = x^{1}$ . This completes the proof of the necessary part. □

The synthesis provided in Theorem III.7 is implicit, as the operator $φ_{u} (T)$ depends on the state trajectory $x (\cdot)$ , and thus on the control input $u$ . This dependence poses practical limitations, even when a solution exists. In applications such as brain stimulation via tDCS, one is often interested in modulating activity over short time horizons. In what follows, we derive an expansion of the matrix $φ_{u} (T)$ in the small-time regime that enables accurate synthesis in small time horizons.

The proof of the following is presented in Section C-B.

Proposition III.8. Under the assumptions of Theorem III.7, the following expansion holds

φ_{u} (T) \underset{T ~ 0}{=} (Id - e^{- T U_{T} (x^{0})}) U_{T} {(x^{0})}^{- 1} D ϕ_{T} (x^{0}) - (e^{T U_{T} (x^{0})} - Id) U_{T} {(x^{0})}^{- 1} + 𝒪 (T^{3})

(23)

In particular, the following expansion also holds

Id - 𝒜_{T} (x^{0}) φ_{u} (T) \underset{T ~ 0}{=} D ϕ_{T} {(x^{0})}^{- 1} e^{- T U_{T} (x^{0})} + 𝒪 (T^{2})

(24)

Remark III.9. It is worth noting that the $𝒪 (T^{3})$ term in (23) becomes $𝒪 (T^{2})$ in (24) due to the relation

𝒜_{T} (x^{0}) 𝒪 (T^{3}) \underset{T ~ 0}{=} 𝒪 (T^{2}),

since $𝒜_{T} (x^{0})$ includes the factor ${(Id - e^{- T U_{T} (x^{0})})}^{- 1}$ , which behaves like $𝒪 (T^{- 1})$ as $T \to 0$ , given that $U_{T} (x^{0})$ is uniformly bounded with respect to $T > 0$ and $x^{0}$ .

The explicit constant control synthesis of this section is summarized in the following theorem, which includes an estimate of the endpoint error.

Theorem III.10. Let $T > 0$ and $(x^{0}, x^{1}) \in {(R^{d})}^{2}$ . Assume that for any $ℓ \in Z, i \frac{2 π ℓ}{T} \notin σ (D N (ϕ_{T} (x^{0})))$ . Then, the solution $x (\cdot)$ to (4) satisfies $x (T) = x^{1}$ , if and only if, $u \in R^{k}$ solves

B u \underset{T ~ 0}{=} {[e^{T U_{T} (x^{0})} - Id]}^{- 1} U_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0})) + 𝒪 (T^{2})

(25)

where $U_{T} (x^{0}) ≔ D N (ϕ_{T} (x^{0}))$ . In particular, let $\tilde{x} (\cdot)$ denote the solution of (4) corresponding to $\tilde{u} \in R^{k}$ satisfying

B \tilde{u} = {[e^{T U_{T} (x^{0})} - Id]}^{- 1} U_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0})) .

(26)

Then, the following endpoint error estimate holds

|\tilde{x} (T) - x^{1}| \underset{T ~ 0}{=} 𝒪 (T^{2}) .

Proof. First, (25) is an immediate consequence of (17), (18) and (24). Next, using the expansions (14), (2), (23), and the definition of $\tilde{u}$ from (26), we obtain

{\tilde{x} (T) = ϕ}_{T} (x^{0}) + \sum_{n = 0}^{\infty} (- 1)^{n} \frac{T^{n + 1} U_{T} {(x^{0})}^{n}}{(n + 1)!} D ϕ_{T} (x^{0}) B \tilde{u} - φ_{u} (T) B \tilde{u} \underset{T ~ 0}{=} [Id - e^{- T U_{T} (x^{0})}] U_{T} {(x^{0})}^{- 1} D ϕ_{T} (x^{0}) B \tilde{u} - [Id - e^{- T U_{T} (x^{0})}] U_{T} {(x^{0})}^{- 1} D ϕ_{T} (x^{0}) B \tilde{u} - (e^{T U_{T} (x^{0})} - Id) U_{T} {(x^{0})}^{- 1} B \tilde{u} + ϕ_{T} (x^{0}) + 𝒪 (T^{2}) \underset{T ~ 0}{=} x^{1} + 𝒪 (T^{2}) .

(27)

The same observation as in Remark III.9 justifies why the $𝒪 (T^{3})$ term in (23) reduces to $𝒪 (T^{2})$ in (27). □

The following result is immediate. It provides a necessary and sufficient condition for (26) to admit at least one solution.

Corollary III.11. Under the hypotheses of Theorem III.7, (26) admits at least one solution $\tilde{u} \in R^{k}$ if and only if

{[e^{T U_{T} (x^{0})} - Id]}^{- 1} U_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0})) \in Im B .

In this case, the least-norm constant control is given by

\tilde{u} = B^{†} {[e^{T U_{T} (x^{0})} - Id]}^{- 1} U_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0})) .

While the previous synthesis applies in a short time, some practical scenarios require reaching the target state over a longer time. In the following, we show how a step function can be constructed from a short-time synthesis to achieve state transfer in a long-time horizon. As noted in Remark II.5, the solution to (4) can still be represented by (9).

We begin with the following implicit synthesis for step-function controls. Its proof, outlined in Section C-C, closely follows the arguments in Proposition III.6 and Theorem III.7.

Theorem III.12. Let $(x^{0}, x^{1}) \in {(R^{d})}^{2}$ and $T \geq τ$ , where $τ > 0$ be such that for any $ℓ \in Z, i \frac{2 π ℓ}{τ} \notin σ (D N (ϕ_{T} (x^{0})))$ . Define the step-function $u_{s f, τ} : [0, T] \to R^{k}$ ,

u_{s f, τ} (t) = \{\begin{array}{l} 0 & if 0 \leq t \leq T - τ \\ u & if T - τ < t \leq T \end{array}

where $u \in R^{k}$ . Then, the solution $x_{τ} (\cdot)$ to (4) corresponding to $u_{s f, τ} (\cdot)$ satisfies $x_{τ} (T) = x^{1}$ , if and only if, $u \in R^{k}$ solves

[Id - 𝒜_{τ, T} (x^{0}) φ_{u} (τ, T)] B u = 𝒜_{τ, T} (x^{0}) (x^{1} - ϕ_{T} (x^{0})) .

(28)

Here, $U_{T} (x^{0}) = D N (ϕ_{T} (x^{0}))$ and

𝒜_{τ, T} (x^{0}) = D ψ_{τ} (ϕ_{T} (x^{0})) {[Id - e^{- τ U_{T} (x^{0})}]}^{- 1} U_{T} (x^{0}) .

Moreover, $φ_{u} (τ, T) = κ_{u} (τ, T) + η_{u} (τ, T)$ , with $κ_{u} (τ, T)$ and $η_{u} (τ, T)$ defined in (15) and (16), respectively, where the integrals are taken over $[T - τ, T]$ .

Then, one has the following synthesis in a long-time horizon using explicit step-function controls.

Theorem III.13. Let $(x^{0}, x^{1}) \in {(R^{d})}^{2}$ and $T \geq τ$ , where $τ > 0$ be such that for any $ℓ \in Z, i \frac{2 π ℓ}{τ} \notin σ (D N (ϕ_{T} (x^{0})))$ . Let ${\tilde{u}}_{τ} \in R^{k}$ be a solution of

B {\tilde{u}}_{τ} = {[e^{τ U_{T} (x^{0})} - Id]}^{- 1} U_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0}))

where $U_{T} (x^{0}) ≔ D N (ϕ_{T} (x^{0}))$ . Then, the solution ${\tilde{x}}_{τ} (\cdot)$ to (4) corresponding to the step function ${\tilde{u}}_{s f, τ} : [0, T] \to R^{k}$ ,

{\tilde{u}}_{s f, τ} (t) = \{\begin{array}{l} 0 & if 0 \leq t \leq T - τ \\ {\tilde{u}}_{τ} & if T - τ < t \leq T \end{array}

satisfies

|{\tilde{x}}_{τ} (T) - x^{1}| \underset{τ ~ 0}{=} 𝒪 (τ^{2}) .

The proof of Theorem III.13 follows the same arguments as Theorem III.10 and is therefore omitted for brevity.

2) Backward nominal-state synthesis: Building on the solution representation of (4) given in Theorem II.4, we aim to construct step-function controls under the assumption that the terminal condition $x (T) = x^{1}$ is achieved. In contrast to the forward nominal-state synthesis, the implicit control derived here does not provide a direct condition guaranteeing that the target state is reached.

We begin with the key result on which the implicit synthesis is based. The proof is presented in Section C-A.

Proposition III.14. Let $T > 0$ and $(x_{0}, u) \in R^{d} \times R^{k}$ . The solution $x (\cdot)$ of (4) can be expanded as

x (t) = ϕ_{t} (x^{0} + \sum_{n = 0}^{\infty} \frac{t^{n + 1} Z_{t}^{n}}{(n + 1)!} P_{t} B u - ζ_{u} (t) B u)

(29)

for all $t \in [0, T]$ . Here $ζ_{u} (t) = ξ_{u} (t) + χ_{u} (t)$ ,

ξ_{u} (t) = \sum_{n = 1}^{\infty} \int_{0}^{t} \frac{s^{n + 1}}{(n + 1)!} [\frac{d}{d s} Z_{s}^{n}] P_{s} d s,

(30)

χ_{u} (t) = \sum_{n = 1}^{\infty} \int_{0}^{t} \frac{s^{n}}{n!} Z_{s}^{n - 1} D^{2} ψ_{s} (x (s)) \dot{x} (s) d s,

(31)

where $P_{t} = D ψ_{t} (x (t))$ , $Z_{t} = D N (ψ_{t} (x (t)))$ and $D^{2} ψ_{s} (x (s))$ is the second derivative of $ψ_{s}$ at $x (s)$ .

Theorem III.15. Let $T > 0$ and $(x^{0}, x^{1}) \in {(R^{d})}^{2}$ . Assume that for any $ℓ \in Z, i \frac{2 π ℓ}{T} \notin σ (D N (ψ_{T} (x^{1})))$ . Any constant control $u \in R^{k}$ that ensures the corresponding solution $x (\cdot)$ of (4) satisfies $x (T) = x^{1}$ necessarily solve

(Id - ℬ_{T} (x^{1}) ζ_{u} (T)) B u = ℬ_{T} (x^{1}) (ψ_{T} (x^{1}) - x^{0}) .

(32)

Here $ζ_{u} (T) \in ℳ_{d} (R)$ is defined in Proposition III.14, and letting $V_{T} (x^{1}) = D N (ψ_{T} (x^{1}))$ , one has

ℬ_{T} (x^{1}) = D ϕ_{T} (ψ_{T} (x^{1})) {(e^{T V_{T} (x^{1})} - Id)}^{- 1} V_{T} (x^{1}) .

(33)

Remark III.16. The implicit backward nominal-state synthesis in Theorem III.15 provides only a necessary condition for a constant control $u \in R^{k}$ to achieve $x (T) = x^{1}$ in (4). In contrast, the forward synthesis (17) yields a sufficient condition: any $u$ solving it guarantees $x (T) = x^{1}$ , and thus also satisfies the backward condition.

Let us now present the proof of Theorem III.15.

Proof of Theorem III.15. Letting $t = T$ in (29), one gets from $x (T) = x^{1}$ and left multiplication by $V_{T} (x^{1})$ that

V_{T} (x^{1}) (ψ_{T} (x^{1}) - x^{0}) = - V_{T} (x^{1}) ζ_{u} (T) B u + (e^{T V_{T} (x^{1})} - Id) D ψ_{T} (x^{1}) B u

(34)

by (2). Since ${[D ψ_{T} (x^{1})]}^{- 1} = D ϕ_{T} (ψ_{T} (x^{1}))$ , left multiplying (34) by ${(e^{T V_{T} (x^{1})} - Id)}^{- 1}$ and $D ϕ_{T} (ψ_{T} (x^{1}))$ yields

D ϕ_{T} (ψ_{T} (x^{1})) {(e^{T V_{T} (x^{1})} - Id)}^{- 1} V_{T} (x^{1}) (ψ_{T} (x^{1}) - x^{0}) = B u - D ϕ_{T} (ψ_{T} (x^{1})) {(e^{T V_{T} (x^{1})} - Id)}^{- 1} V_{T} (x^{1}) ζ_{u} (T) B u .

It follows that $u \in R^{k}$ solves the implicit equation

[Id - ℬ_{T} (x^{1}) ζ_{u} (T)] B u = ℬ_{T} (x^{1}) (ψ_{T} (x^{1}) - x^{0})

which is (32), where $ℬ_{T} (x^{1})$ is defined by (33). □

Using the same machinery that leads to Theorem III.10, we obtain the following result. The proof is omitted for brevity.

Theorem III.17. Under the assumptions of Theorem III.15, let $\overline{x} (\cdot)$ denote the solution of (4) corresponding to a constant control $\overline{u} \in R^{k}$ satisfying

B \overline{u} = {[Id - e^{- T V_{T} (x^{1})}]}^{- 1} V_{T} (x^{1}) (ψ_{T} (x^{1}) - x^{0})

(35)

where $V_{T} (x^{1}) ≔ D N (ψ_{T} (x^{1}))$ . Then, the following endpoint error estimate holds

|\overline{x} (T) - x^{1}| \underset{T ~ 0}{=} 𝒪 (T^{2}) .

The following result provides a necessary and sufficient condition for (35) to admit at least one solution.

Corollary III.18. Under the hypotheses of Theorem III.15, (35) admits at least one solution $\tilde{u} \in R^{k}$ if and only if

{[Id - e^{- T V_{T} (x^{1})}]}^{- 1} V_{T} (x^{1}) (ψ_{T} (x^{1}) - x^{0}) \in Im B .

In this case, the least-norm constant control is given by

\overline{u} = B^{†} {[Id - e^{- T V_{T} (x^{1})}]}^{- 1} V_{T} (x^{1}) (ψ_{T} (x^{1}) - x^{0}) .

Finally, one has the following result that synthesizes a step function for a large time horizon.

Theorem III.19. Let $(x^{0}, x^{1}) \in {(R^{d})}^{2}$ and $T \geq τ$ , where $τ > 0$ be such that for any $ℓ \in Z, i \frac{2 π ℓ}{τ} \notin σ (D N (ψ_{T} (x^{1})))$ . Let ${\overline{u}}_{τ} \in R^{k}$ be a solution of

B {\overline{u}}_{τ} = {[Id - e^{- τ V_{T} (x^{1})}]}^{- 1} V_{T} (x^{1}) (ψ_{T} (x^{1}) - x^{0})

where $V_{T} (x^{1}) = D N (ψ_{T} (x^{1}))$ . Then, the solution $x (\cdot)$ to (4) corresponding to the step function ${\overline{u}}_{s f} : [0, T] \to R^{k}$ ,

{\overline{u}}_{s f} (t) = \{\begin{array}{l} 0 & if 0 \leq t \leq T - τ \\ \overline{u} & if T - τ < t \leq T \end{array}

satisfies

|x (T) - x^{1}| \underset{τ ~ 0}{=} 𝒪 (τ^{2}) .

Remark III.20. A key distinction between the forward and backward nominal-state syntheses lies in their respective spectral conditions: the forward synthesis depends on the spectrum at the initial state $x^{0}$ , while the backward synthesis relies on that at the target state $x^{1}$ .

In practice, only one of these conditions may be satisfied, depending on the region of the state space. This asymmetry underscores the complementarity of the two syntheses—each offers a valid control strategy under different local spectral properties, making their coexistence practically valuable.

This contrast also reflects a deeper difference between linear and nonlinear control. In the linear case, controllability depends on a uniform spectral condition (cf. Proposition III.2); failure at one point implies failure everywhere. In the nonlinear setting, local spectral conditions at either the initial or target can independently ensure controllability, highlighting a key advantage of the nonlinear synthesis framework.

C. Forward and Backward syntheses

Although Sections III-B1 and III-B2 focus on the forward and backward nominal-state syntheses, respectively, it is worth emphasizing that both arise from a broader family of synthesis strategies derived from the dual representations of the nonlinear system (4) in Theorems II.2 and II.4. Each representation admits at least two distinct expansion strategies—depending on the chosen integration by parts identity—that lead to meaningful control synthesis formulations.

Specifically, starting from the solution representation (7), we observe two synthesis paths:

Applying the identity $\int_{0}^{T} D ψ_{t} d t = \int_{0}^{T} t^{'} D ψ_{t} d t$ yields the implicit backward nominal-state synthesis fully investigated in Section III-B2.
Using instead $\int_{0}^{T} D ψ_{t} d t = \int_{0}^{T} (t - T)^{'} D ψ_{t} d t$ gives rise to what we refer to as an implicit backward initial-state synthesis. This formulation yields an expression of the control that provides a necessary and sufficient condition for reachability, and it is not analyzed in detail here.

Similarly, from the solution representation (9), two analogous options arise:

Using the identity $\int_{0}^{T} D ϕ_{T - t} d t = \int_{0}^{T} (t - T)^{'} D ϕ_{T - t} d t$ leads to the implict forward nominal-state synthesis fully developed in Section III-B1.
Alternatively, using $\int_{0}^{T} D ϕ_{T - t} d t = \int_{0}^{T} t^{'} D ϕ_{T - t} d t$ results in what we refer to as an implicit forward final-state synthesis. This yields a formally valid expansion and a necessary condition on the control, and it is not investigated in this work.

D. Reachability via control synthesis methods

In this section, we summarize the main results from Sections III-B1 and III-B2 into a set of operational insights framed in terms of the input matrix $B$ . Specifically, we analyze how the set of states reachable from a given initial condition $x^{0} \in R^{d}$ over the interval [ $0, T$ ] depends on the structure of $B$ , under both constant and step-function control strategies synthesized via the proposed methods.

Given $x^{0} \in R^{d}$ , we say that a state $x^{1} \in R^{d}$ is c-reachable over $[0, T]$ from $x^{0}$ if there exists a constant control $u \in R^{k}$ such that the solution $x (\cdot)$ of (4) corresponding to $u$ satisfies $x (T) = x^{1}$ . The c-reachable set is denoted by (for step-function controls, we denote it by $ℛ_{s c} (τ, T, x^{0})$ )

ℛ_{c} (T, x^{0}) = \{x (T) |\begin{array}{l} x (\cdot) solves (4) \\ with constant control u \in R^{k} \end{array}\} .

(36)

We say that (4) is controllable over $[0, T]$ from $x^{0}$ with a constant control if $ℛ_{c} (T, x^{0}) = R^{d}$ . If this holds for every $x^{0} \in R^{d}$ , then (4) is said to be completely controllable over $[0, T]$ with constant controls.

1) The linear case: It is convenient to begin with the linear system (10), which helps understanding the fully nonlinear setting. In the linear control framework, controllability with a time-varying control is guaranteed if the Kalman rank condition or, equivalently, the invertibility of the controllability Gramian holds. These conditions are independent of $T > 0$ and $x^{0}$ , and they allow for flexibility in choosing the input matrix $B \in ℳ_{d, k} (R)$ satisfying what we refer to as the time-varying controllability condition.

However, as shown in Remark III.1, the question becomes more delicate when restricting to constant controls. Notably, even the trivial case $B = Id$ may fail to achieve controllability with a constant control, despite satisfying the time-varying controllability condition.

Under the spectral condition in Proposition III.2, Corollary III.3 shows that (10) is completely controllable with constant controls over $[0, T]$ if and only if $k = d$ and $B$ is invertible. Otherwise, if $k < d$ or $B$ is not invertible, then

ℛ_{c} (T, x^{0}) ⊈ R^{d}, for all x^{0} \in R^{d} .

Hence, the size of $ℛ_{c} (T, x^{0})$ is determined by the rank of : the smaller $k$ is relative to $d$ , the fewer the number of target states that are reachable from $x^{0}$ via constant control.

2) The Nonlinear Case: We now turn to the fully nonlinear setting and analyze the c-reachable set $ℛ_{c} (T, x^{0})$ , leveraging the forward nominal-state synthesis from Section III-B1.

Unlike the linear case, the synthesis equation in Theorem III.7 is implicit in $u$ unless $φ_{u} (T) = 0$ . When $φ_{u} (T) = 0$ , the analysis from Section III-D1 directly applies. If $φ_{u} (T) \neq$ 0, Proposition III.8 shows that for sufficiently small $T > 0$ , the control $u$ satisfies (25). As a result, (4) is completely controllable with constant controls over $[0, T]$ for small $T > 0$ if and only if $k = d$ and $B$ is invertible. Otherwise, when $k < d$ or $B$ is not invertible, the c-reachable set satisfies

ℛ_{c} (T, x^{0}) ⊈ R^{d}, for all x^{0} \in R^{d},

and its size is again governed by the rank of $B$ .

Still, when $φ_{u} (T) \neq 0$ , the small-time reachability result allows one to extend the synthesis to arbitrary horizons $T > 0$ via step-function controls. Specifically, let $τ > 0$ be such that for any $ℓ \in Z, \frac{2 π ℓ}{τ} \notin σ (D N (ϕ_{T} (x^{0}))$ . Then, for $T \leq τ$ , we have $ℛ_{s c} (τ, T, x^{0}) = ℛ_{c} (τ, x^{0})$ . While for $T > τ$ , Theorem III.12 shows that system (4) is completely controllable with step-function controls over $[0, T]$ if and only if $k = d$ and $B$ is invertible. Otherwise, one has

ℛ_{s c} (τ, T, x^{0}) ⫋ R^{d}, for all x^{0} \in R^{d}

and the sc-reachable set structure again depends on $rank (B)$ .

Let us characterize $ℛ_{c} (T, x^{0})$ for the input matrix traditionally considered in network neuroscience [16], namely

B = [\begin{array}{l} e_{1} & e_{2} & \dots & e_{k} \end{array}], k < d,

(37)

where $e_{i}$ is the $i$ -th canonical basis vector of $R^{d}$ .

Proposition III.21. Let $T > 0$ and $(x^{0}, x^{1}) \in {(R^{d})}^{2}$ . Assume that for any $ℓ \in Z, i \frac{2 π ℓ}{T} \notin σ (U_{T} (x^{0}))$ . If the input matrix $B$ is given by (37), then the c-reachable set satisfies

ℛ_{c} (T, x^{0}) \underset{T ~ 0}{=} ϕ_{T} (x^{0}) + ker M_{T} (x^{0}) + 𝒪 (T^{2}), dim (ker M_{T} (x^{0})) = k,

(38)

and the constant input $u \in R^{k}$ that steers (4) from $x^{0} \in R^{d}$ to $x^{1} \in ϕ_{T} (x^{0}) + ker M_{T} (x^{0})$ is given by

u = V_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0})) .

(39)

Here $U_{T} (x^{0}) ≔ D N (ϕ_{T} (x^{0})), V_{T} (x^{0}) \in ℳ_{k, d} (R)$ collects the first $k$ rows of $ℬ_{T} (x^{0}) ≔ {[e^{T U_{T} (x^{0})} - Id]}^{- 1} U_{T} (x^{0})$ , and $M_{T} (x^{0}) \in ℳ_{d - k, d} (R)$ collects its remaining $d - k$ rows.

Proof. First, it follows from Theorem III.10 that the solution $x (\cdot)$ to (4) corresponding to a constant control $u \in R^{k}$ solving

B u = {[e^{T U_{T} (x^{0})} - Id]}^{- 1} U_{T} (x^{0}) (x^{1} - ϕ_{T} (x^{0}))

(40)

satisfies (see (27))

x (T) \underset{T ~ 0}{=} x^{1} + 𝒪 (T^{2}) .

(41)

Next, with $B$ as in (37), one deduces from (40) that

[\begin{matrix} {Id}_{k} \\ O_{d - k, k} \end{matrix}] u = [\begin{matrix} V_{T} (x^{0}) \\ M_{T} (x^{0}) \end{matrix}] (x^{1} - ϕ_{T} (x^{0}))

(42)

where ${Id}_{k} \in ℳ_{k} (R)$ is the identity matrix, and $O_{d - k, k} \in ℳ_{d - k, k} (R)$ is the zero matrix. Since $ℬ_{T} (x^{0})$ is invertible, one has $rank M_{T} (x^{0}) = d - k$ and therefore $dim (ker M_{T} (x^{0})) = d - (d - k) = k$ . Finally, (42) is equivalent to (39) and

x^{1} \in ϕ_{T} (x^{0}) + ker M_{T} (x^{0}) .

(43)

Combining (36), (41) and (43) yields (38). □

Remark III.22. In Proposition III.21, the first-order characterization of the c-reachable set as an affine subspace is valid as long as the desired target shift $x^{1} - ϕ_{T} (x^{0}) = 𝒪 (T)$ . This ensures that the required constant input (39) remains bounded, since to leading order $T \to 0$ , one has $| u | \approx |x^{1} - ϕ_{T} (x^{0})| / T$ .

Using the thin QR factorization of $M_{T} {(x^{0})}^{⊤} \in R^{d \times (d - k)}$ ,

M_{T} {(x^{0})}^{⊤} = [\begin{array}{l} Q_{1} & Q_{2} \end{array}] [\begin{matrix} R_{1} \\ 0 \end{matrix}], \{\begin{array}{l} Q_{1} \in R^{d \times (d - k)}, Q_{2} \in R^{d \times k}, \\ R_{1} \in R^{(d - k) \times (d - k)}, \end{array}

one finds that $ker M_{T} (x^{0}) = span (Q_{2})$ .

Remark III.23. First, the c-reachable set characterization and the associated constant input synthesis in Proposition III.21 for a small time horizon $T > 0$ when the input matrix $B$ is given by (37) can be improved by explicitly considering the second derivative $D^{2} ϕ_{T - t}$ in the expansion. We defer this analysis to future works. Next, let $T > 0$ and $x^{0} \in R^{d}$ be such that $D N (ϕ_{T} (x^{0}))$ satisfies the spectral condition in Theorem III.7. Future work should investigate the structure and effective dimension of the c-reachable sets $ℛ_{c} (T, x^{0})$ as functions of $rank (B)$ for a general input matrix $B \in ℳ_{k, d} (R)$ .

IV. Comments on the main results of Section III

First, if the spectral norm condition $‖ W ‖ < λ_{min} (D)$ is satisfied, then the eigenvalue-related assumptions required in Proposition III.2, Theorems III.7, and III. 15 are automatically fulfilled, as discussed, e.g., in Remark III.4. Note that this condition ensures that the nonlinear system (4) is contracting.

Contracting dynamics possess a range of desirable properties–see [20]–and contraction has become a standard structural assumption in many recent works on recurrent neural networks, such as [7], [11]. While the condition on $‖ W ‖$ is a sufficient criterion for contraction, it serves here as a simple and practically verifiable condition that ensures well-posedness of the synthesis proposed in this work.

More generally, the control syntheses developed in this paper apply to any nonlinear system of the form $\dot{x} (t) = N (x (t)) + B u$ , where $u \in R^{k}$ is constant and $N \in C^{2} (R^{d}; R^{d})$ satisfies the regularity conditions

sup_{x \in R^{d}} ‖ D N (x) ‖ \leq C, sup_{x \in R^{d}} ‖D^{2} N (x)‖ \leq C

for some constant $C > 0$ . Notably, we do not require any explicit boundedness condition on the vector field $N$ itself, making the framework broadly applicable.

This flexibility is particularly relevant for applications in machine learning, where recurrent neural networks (RNNs) often employ unbounded and non-smooth activation functions such as the $ReLU : f_{i} (x_{i}) = ReLU (x_{i}) = max (0, x_{i})$ . ReLU is differentiable everywhere except at $x_{i} = 0$ , and it is globally Lipschitz on $R$ . In this case, the solution representations (7) and (9) no longer apply, and our control syntheses cannot be directly used. To overcome this limitation, one can consider smooth approximations of the ReLU function (e.g., softplus [14], swish function [29]), which restore differentiability and allow the application of our synthesis results.

V. Examples and numerical simulations

This section presents some numerical simulations that underpin our theoretical study. The pipeline of the experiments and analysis is visualized in Fig. 2. We included three major types of models: linear systems, “vanilla” tanh RNNs, and MINDy models [32]. While vanilla RNNs are widely used in theoretical neuroscience as task-performing models, MINDy models represent another branch of models that directly approximate experimental neural data, demonstrating the different applicative contexts of our theoretical framework.

Fig. 2. — Pipeline of the numerical experiments.

A. Implementation

We implement the linear system control synthesis (13), the nonlinear forward nominal-state synthesis (26), and the nonlinear backward nominal-state synthesis (35) in python. The scripts will be available after publication. For the nonlinear synthesis, the flow is numerically integrated using the odeint function from the torchdiffeq package [9]. We adopted the Runge-Kutta method of order 7(8) of Dormand-Prince-Shampine, which is the highest-order method available from the package. The method accepts a new step if $RMS (ε) < 10^{- 14} + 10^{- 13} * RMS (x)$ , where $RMS (\cdot)$ represents the root mean square norm, $ε$ being the estimated error and $x$ the current state. We found that this setting provided reasonable control of the numerical error related to flow integration. The Jacobian matrices were computed through auto-differentiation in pytorch. Experiments were conducted on a desktop computer with an Nvidia RTX 3080 GPU. The run time for each synthesis increased as the time horizon $T$ and number of features (neurons) $d$ increased, but remained within the scope of several hundred milliseconds to several seconds. For a fixed $T$ and multiple desired trajectories (i.e., different pairs of $(x^{0}, x^{1})$ ), the synthesis is efficiently parallelized, and the run time remains almost constant when increasing the number of trajectories from one to 1000.

To construct a $d$ -dimensional linear system, we randomly sampled the entries of $W$ from the normal distribution $𝒩 (0,1 / d)$ . The decay matrix was set to $(λ_{max} (W) - λ_{0}) Id$ , where $λ_{max} (W) \in R$ is the maximum of the real parts of the eigenvalues of $W$ and $λ_{0} \in R$ is the desired maximum real part of the eigenvalues of the dynamic matrix $A = W - D$ . We tried both $λ_{0} = - 0.1$ and $λ_{0} = 0.1$ to construct stable and unstable systems respectively.

To construct a “vanilla” RNN, we set $D = Id$ and the activation function to tanh. We followed [35] and imposed a “random plus low rank” structure on the connectivity matrix $W = J + {m n}^{⊤}$ . The matrix $J \in ℳ_{d} (R)$ was sampled from the normal distribution $𝒩 (0, g^{2} / d)$ , where $g$ is a scaling factor. $m, n \in ℳ_{d, p} (R)$ determines the rank of the low-rank component. Such a type of model appears frequently in theoretical neuroscience studies [23]. We considered three subtypes of vanilla RNNs. The first subtype is referred to as ‘small norm tanh RNN’, constructed with $g = 0.5$ and $p = 0$ (no low-rank component). We found that such RNNs satisfied $‖ W D f (0) ‖ \approx 1 = λ_{min} (D)$ . For the next two subtypes, $g$ was set to 0.9 [35] and the spectral norm condition was violated. $m \in R^{d}$ was sampled from standard normal distribution. $n \in R^{d}$ was sampled from normal distribution $𝒩 (0,1 / d^{2})$ for the second subtype and set as $n = \frac{1.1}{d} m$ for the third subtype. It can be shown that systems of the second subtype will be monostable while systems of the third subtype will be bistable [35] as the number of neurons $d$ tends to infinity.

To further evaluate the method on realistic models fit to experimental data, we also included a set of Mesoscale Individualized Neurodynamic (MINDy) models from [32]. MINDy models contained 100 interconnected units representing 100 brain areas, and the parameters were optimized to approximate the activation time series of these areas measured through functional magnetic resonance imaging (fMRI). Unlike most RNNs, the activation function of MINDy is heterogeneous: $f (s) = \sqrt{α^{2} + (b s + 0.5)^{2}} - \sqrt{α^{2} + (b s - 0.5)^{2}}$ , where $b = 20 / 3$ is fixed and $α$ is optimized over the data and differs across the 100 units. It is also worth noting that the origin is unstable in most of the MINDy models, indicating that the spectral norm condition $‖ W D f (0) ‖ < λ_{min} (D)$ was not met. However, we found that the models satisfy the eigenvalue condition, which enabled the synthesis.

In the first experiment, we analyzed the endpoint error of the controlled trajectory for $T \in \{2^{- 2}, 2^{- 1}, \dots, 2^{6}\}$ when $B = Id$ . For each $T$ , we randomly generated 5 stable and 5 unstable linear systems, 5 small-norm, 5 monostable and 5 bistable tanh RNNs, and randomly selected 5 fitted MINDy models (from a pool of 106 models). All models were 100-dimensional. Then, we randomly generated 40 initial states $x^{0} ~ 𝒩 (0, Id)$ . To investigate how the deviation of $x^{1}$ from the autonomous flow end point $ϕ_{T} (x^{0})$ influences the performance of the method, we set $x^{1} = ϕ_{T} (x^{0}) + ε$ where $ε ~ 𝒩 (0, σ^{2} \cdot Id$ ). In half of the trials, $σ^{2}$ was set to 0.1 (“small deviation”); in the other half, $σ^{2}$ was set to 0.5 (“large deviation”). Control input was computed using the nonlinear forward (26) and backward (35) syntheses, and the linear synthesis (13). We recall that for the nonlinear system (4), when $B = Id$ , the linearized system at $x^{0}$ is given by

\dot{x} (t) = D N (x^{0}) x (t) + N (x^{0}) + u .

Using (13), the following input $u$ drives the linearized system from $x^{0}$ to $x^{1}$ in exact time $T$

u = {(e^{T A} - Id)}^{- 1} A (x^{1} - e^{T A} x^{0}) - N (x^{0})

(44)

where $A = D N (x^{0})$ . We have also tried to linearize the systems at the origin (which is always a fixed point) and obtained qualitatively similar results.

In the second experiment, we followed the tradition in network neuroscience [16] and considered the input matrix defined by (37). For simplicity, we focused on 128-dimensional “small norm” tanh RNNs (with $g = 0.5$ and $p = 0$ ) and $T \in {0.25,0.5,1, 2,4}$ . We considered $k \in {1,64,96,120,126,128}$ , ranging from an extremely underactuated system to a fully actuated system. We randomly generated 5 RNNs and selected 20 initial states $x^{0} ~ 𝒩 (0, Id)$ for each model and each combination of $T$ and $k$ . Here, $x^{1} = ϕ_{T} (x^{0}) + Q_{2} ξ$ by Proposition III.21, and $Q_{2}$ following Remark III.22, with entries of $ξ$ drawn from $𝒩 (0,0.01)$ .

B. Results

Results of the first experiment were summarized in Fig. 3. Monostable and bistable tanh RNNs were combined into “Big norm tanh RNN” due to similarity. For linear systems, the error remained small. For nonlinear systems, error generally increased as the horizon $T$ increased. The forward nominal state synthesis performed better than the backward nominal state synthesis, and both of them generally performed better than linearization, particularly for large $T$ . The differences were more evident when $x^{1} - ϕ_{T} (x^{0})$ is small (with standard deviation 0.1, left panels), where the forward synthesis sometimes outperformed linearization by orders of magnitude. Note that the MINDy models were trained to approximate timeseries with unit variance, so a “small deviation” with standard deviation 0.1 already represents a physically sizable difference, indicating the practical significance of our method.

Results of the second experiment were summarized in Fig. 4. The nonlinear synthesis was still much better than linearization in underactuated systems, as expected. However, surprisingly, the endpoint error did not change dramatically and even decreased as the number of actuators decreased, as long as $x^{1}$ is close to the reachable set. These findings suggest that Proposition III.21 did provide a good estimation of the reachable set, and that the synthesis worked in underactuated systems just as good as (if not better than) fully-actuated systems as long as $x^{1}$ is reachable.

VI. Concluding remarks and perspectives

This paper addressed the control synthesis problem for a class of nonlinear Hopfield-type recurrent neural networks motivated by neurostimulation applications. Using a solution representation that generalizes the variation of constants formula, we derived constant and piecewise constant inputs capable of steering the network to a desired target within a prescribed time interval. For linear activation functions, the exact controllability with constant input is possible over arbitrary time horizons. For a small time horizon, we also provided a characterization of the reachable set when the input matrix $B$ directly actuates a subset of nodes, showing that its dimension equals the rank of $B$ , and that its basis can be computed efficiently via a thin QR factorization. Moreover, the constant input that guarantees reachability is given by a closed-form algebraic condition, which makes the synthesis directly applicable for tDCS.

Future work will explore the synthesis of time-varying control inputs based on the proposed framework, with an emphasis on robustness to parameter uncertainty and external disturbances. As large-scale systems like (4) are highly sensitive to such perturbations, designing controls that adapt dynamically while minimizing energy costs remains a key challenge—one that links controllability with energetic performance.

Supplementary Material

NIHMS2152725-supplement-1.pdf^{(122.5KB, pdf)}

Acknowledgments

This work is partially supported by grant R21MH132240 from the US National Institutes of Health to SC.

Biographies

Cyprien Tamekue was born in Nkongsamba, Cameroon, in November 1994. He received his B.Sc. degree in Fundamental Computing from the University of Douala, Cameroon, in July 2016. He obtained a M.Sc. degree in Fundamental Mathematics and Applications from Institut de Mathématiques et de Sciences Physiques, Benin, in 2019 and a M.Sc. degree in Theoretical and Applied Mathematics from Université Paris Dauphine-PSL in 2020. He completed his Ph.D. in Control Theory and Mathematical Neuroscience at Université Paris-Saclay in October 2023.

Dr. Tamekue is currently a Postdoctoral Research Associate at Washington University in St. Louis, McKelvey School of Engineering, in the Brain Dynamics and Control Research Group. He was previously a Postdoctoral Researcher in Control Theory at Institut Polytechnique des Sciences Avancées, INRIA, L2S, France. Dr. Tamekue’s research focuses on the mathematical modeling and control of nonlinear dynamical systems in large-scale brain networks, neural field models, and PDEs from sub-Riemannian geometry. He has contributed to modeling visual phenomena, such as the visual MacKay effect, by leveraging ideas from control theory to neural field models.

Ruiqi Chen received the B.Sc. degree in Intelligence Science and Technology from Peking University, China, 2021. He is currently working towards a Ph.D. degree in Neurosciences at Washington University in St. Louis, USA, under the supervision of Professor ShiNung Ching and Professor Todd Braver.

His research focuses on the modeling and analysis of brain dynamics supporting high-level cognitive functions, and the development of tools to analyze large-scale dynamical systems with the combination of theory-driven and data-driven approaches.

ShiNung Ching (Senior Member, IEEE), is Professor in the Department of Electrical and Systems Engineering at Washington University in St. Louis (St. Louis, USA). His research interests are at the intersection of engineering and computational neuroscience, particularly in using systems and control theoretic concepts to study the links between dynamics and function in neuronal networks. Dr. Ching completed his B.Eng (Hons.) and M.A.Sc degrees in Electrical and Computer Engineering from McGill University, Canada and the University of Toronto, Canada. He earned his Ph.D. in Electrical Engineering from the University of Michigan in 2009 in the area of systems and control theory. He subsequently completed post-doctoral training in computational neuroscience and anesthesiology at the Massachusetts Institute of Technology and the Harvard Medical School. He is an author on over 100 publications in academic journals and conferences and the textbook ‘Quasilinear Control’. Dr. Ching has received the CAREER Award from the US National Science Foundation, the Young Investigator Program award from the US AFOSR, a Career Award at the Scientific Interface from the Burroughs-Wellcome Fund, and has served as principal investigator on research projects funded by the National Institutes of Health, the National Science Foundation and the US Department of Defense.

Contributor Information

Cyprien Tamekue, Department of Electrical and Systems Engineering, Washington University in St. Louis, St. Louis, 63130, MO, USA.

Ruiqi Chen, Neurosciences Program in the Division of Biology and Biomedical Sciences, Washington University in St. Louis, St. Louis, 63130, MO, USA.

ShiNung Ching, Department of Electrical and Systems Engineering, Washington University in St. Louis, St. Louis, 63130, MO, USA.

References

[1].Agračev AA, Barilari D, and Boscain U, A comprehensive introduction to sub-Riemannian geometry. Vol 181, Cambridge University Press, 2019. [Google Scholar]
[2].Agračev AA, and Gamkrelidze RV, “ The exponential representation of flows and the chronological calculus,” Mathematics of the USSR-Sbornik, vol. 35, no. 6, pp. 727, 1978. [Google Scholar]
[3].Albertini F and Sontag ED, “For neural networks, function determines form,” Neural Networks, vol. 6, no. 7, pp. 975–990, 1993. [Google Scholar]
[4].Bergmann TO, Groppa S, Seeger M, Mölle M, Marshall L, and Siebner HR, “Acute changes in motor cortical excitability during slow oscillatory and constant anodal transcranial direct current stimulation,” Journal of Neurophysiology, vol. 102, no. 4, pp. 2303–2311, 2009. [DOI] [PubMed] [Google Scholar]
[5].Breakspear M, “ Dynamic models of large-scale brain activity,” Nature Neuroscience, vol. 20, no. 3, pp. 340–352, 2017. [DOI] [PubMed] [Google Scholar]
[6].Bressan A and Piccoli B, Introduction to the Mathematical Theory of Control, Vol. 1, American Institute of Mathematical Sciences, Springfield, 2007. [Google Scholar]
[7].Centorrino V, Gokhale A, Davydov A, Russo G, and Bullo F, “Euclidean contractivity of neural networks with symmetric weights,” IEEE Control Systems Letters, vol. 7, pp. 1724–1729, 2023. [Google Scholar]
[8].Chen R, Singh M, Braver TS, and Ching S, “Dynamical Models Reveal Anatomically Reliable Attractor Landscapes Embedded in Resting State Brain Networks,” Imaging Neuroscience, 2025, doi: 10.1162/imag_a_00442. Available: https://doi.org/10.1162/imag_a_00442. [DOI] [Google Scholar]
[9].Chen RTQ, Rubanova Y, Bettencourt J, and Duvenaud D, “Neural Ordinary Differential Equations,” Advances in Neural Information Processing Systems. 2018. [Google Scholar]
[10].Coron JM, Control and nonlinearity. No. 136, American Mathematical Soc, 2007. [Google Scholar]
[11].Davydov A, Proskurnikov AV and Bullo F (2022, June). “Non-Euclidean contractivity of recurrent neural networks,” In 2022 American Control Conference (ACC) (pp. 1527–1534). IEEE. [Google Scholar]
[12].DeGiorgio CM, and Krahl SE, “Neurostimulation for drug-resistant epilepsy,” Continuum: Lifelong Learning in Neurology, vol. 19, no. 3, pp. 743–755, 2013. [DOI] [PMC free article] [PubMed] [Google Scholar]
[13].Forti M and Tesi A, “New conditions for global stability of neural networks with application to linear and quadratic programming problems,” IEEE Transactions on Circuits and Systems I: Fundamental theory and applications, vol. 47, no. 7, pp. 354–366, 1995. [Google Scholar]
[14].Glorot X, Bordes A, and Bengio Y, “Deep sparse rectifier neural networks,” in Proc. 14th Int. Conf. on Artificial Intelligence and Statistics (AISTATS), pp. 315–323, 2011. [Google Scholar]
[15].Grover S, Wen W, Viswanathan V, Gill CT, and Reinhart RMG, “Long-lasting, dissociable improvements in working memory and long-term memory in older adults with repetitive neuromodulation,” Nature Neuroscience, vol. 25, no. 9, pp. 1237–1246, 2022. [DOI] [PMC free article] [PubMed] [Google Scholar]
[16].Gu S, Pasqualetti F, Cieslak M, Telesford QK, Yu AB, Kahn AE, Medaglia JD, Vettel JM, Miller MB, and Grafton ST, and et al. , “Controllability of structural brain networks,” Nature Communications, vol. 6, no. 1, pp. 8414, 2015. [Google Scholar]
[17].Hirsch MW and Smale S, Differential Equations, Dynamical Systems, and Linear Algebra, Academic Press, New York-London, 1974. Pure and Applied Mathematics, Vol. 60. [Google Scholar]
[18].Hopfield JJ, “Neurons with graded response have collective computational properties like those of two-state neurons,” Proceedings of the National Academy of Sciences, vol. 81, no. 10, pp. 3088–3092, 1984. [Google Scholar]
[19].Liu Y-Y, Slotine J-J, and Barabási A-L, “Controllability of complex networks,” Nature, vol. 473, no. 7346, pp. 167–173, 2011. [DOI] [PubMed] [Google Scholar]
[20].Lohmiller W and Slotine JJE. “On contraction analysis for nonlinear systems,” Automatica, vol. 34 no. 6, pp 683–696, 1998. [Google Scholar]
[21].Kawski M, Chronological calculus in systems and control theory, Mathematics of Complexity and Dynamical Systems, pp. 88, 2011. [Google Scholar]
[22].López-Alonso V, Cheeran B, Río-Rodríguez D, and Fernández-del-Olmo M, “Inter-individual variability in response to non-invasive brain stimulation paradigms”, Neuron, vol. 7, no. 3, pp. 372–380, 2014. [Google Scholar]
[23].Mastrogiuseppe F, and Ostojic S, “Linking connectivity, dynamics, and computations in low-rank recurrent neural networks”, Neuron, vol. 99, no. 3, pp. 609–623, 2018. [DOI] [PubMed] [Google Scholar]
[24].Medaglia JD, “Clarifying cognitive control and the controllable connectome,” Wiley Interdisciplinary Reviews: Cognitive Science, vol. 10, no. 1, p. e1471, 2019. [DOI] [PMC free article] [PubMed] [Google Scholar]
[25].Muldoon SF, Pasqualetti F, Gu S, Cieslak M, Grafton ST, Vettel JM, and Bassett DS, “Stimulation-based control of dynamic brain networks”, PLoS Computational Biology, vol. 12, no. 9, pp. e1005076, 2016. [DOI] [PMC free article] [PubMed] [Google Scholar]
[26].Power JD, Schlaggar BL, Lessov-Schlaggar CN, and Petersen SE, “Evidence for hubs in human functional brain networks”, Neuron, vol. 79, no. 4, pp. 798–813, 2013. [DOI] [PMC free article] [PubMed] [Google Scholar]
[27].Qin Y, Bassett DS, and Pasqualetti F, “Vibrational control of cluster synchronization: Connections with deep brain stimulation,” in Proc. 61st IEEE Conf. on Decision and Control (CDC), pp. 655–661, 2022. [Google Scholar]
[28].Qin Y, Nobili AM, Bassett DS, and Pasqualetti F, “Vibrational stabilization of cluster synchronization in oscillator networks,” IEEE Open Journal of Control Systems, vol. 2, pp. 439–453, 2023. [Google Scholar]
[29].Ramachandran P, Zoph B, and Le QV, “Searching for activation functions,” arXiv preprint, arXiv:1710.05941, 2017. [Google Scholar]
[30].Sarychev AV, “Lie extensions of nonlinear control systems” Journal of Mathematical Sciences, vol. 135, no. 4, pp. 3195–3223, 2006. [Google Scholar]
[31].Singh MF, Wang A, Cole M, Ching S, Braver TS, “Enhancing task fMRI preprocessing via individualized model-based filtering of intrinsic activity dynamics”, NeuroImage, vol. 247, pp. 118836, 2022. [DOI] [PMC free article] [PubMed] [Google Scholar]
[32].Singh MF, Braver TS, Cole MW, and Ching S, “Estimation and validation of individualized dynamic brain models with resting-state fMRI”, NeuroImage, vol. 221, pp. 117046, 2020. [DOI] [PMC free article] [PubMed] [Google Scholar]
[33].Singh MF, Matthew, Wang A, Anxu, Braver TS, and Ching S, “Scalable surrogate deconvolution for identification of partially-observable systems and brain modeling”, Journal of neural engineering, vol. 17, no. 4, pp. 046025, 2020. [DOI] [PMC free article] [PubMed] [Google Scholar]
[34].Schiff SJ, Neural control engineering: the emerging intersection between control theory and neuroscience, MIT Press, 2011. [Google Scholar]
[35].Schuessler F, Dubreuil A, Mastrogiuseppe F, Ostojic S, and Barak O, “Dynamics of random recurrent networks with correlated low-rank structure”, Physical Review Research, vol. 2, no. 1, pp. 013111, 2020. [Google Scholar]
[36].Sontag ED, and Qiao Y, “Remarks on controllability of recurrent neural networks”, Proceedings of the IEEE Conference on Decision and Control., Vol. 1. Institute of Electrical and Electronics Engineers Inc., 1998. [Google Scholar]
[37].Tamekue C and Ching S, “Control analysis and synthesis for general control-affine systems,” arXiv preprint, arXiv:2503.05606, 2025. [Google Scholar]
[38].Tamekue C, Prandi D, and Chitour Y, “A mathematical model of the visual MacKay effect”, SIAM Journal on Applied Dynamical Systems, vol. 23, no. 3, pp. 2138–2178, 2024. [Google Scholar]
[39].Tucsnak M and Weiss G, Observation and control for operator semigroups, Springer, 2009. [Google Scholar]
[40].Yi Z, Heng P-A, and Fu AW-C, “Estimate of exponential convergence rate and exponential stability for neural networks”, IEEE Transactions on Neural Networks, vol. 10, no. 6, pp. 1487–1493, 1999. [DOI] [PubMed] [Google Scholar]
[41].Zhang H, Wang Z, and Liu D, “A comprehensive review of stability analysis of continuous-time recurrent neural networks,” IEEE Transactions on Neural Networks and Learning Systems, vol. 25, no. 7, pp. 1229–1262, 2014. doi: 10.1109/TNNLS.2014.2317880. [DOI] [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

NIHMS2152725-supplement-1.pdf^{(122.5KB, pdf)}

[R1] [1].Agračev AA, Barilari D, and Boscain U, A comprehensive introduction to sub-Riemannian geometry. Vol 181, Cambridge University Press, 2019. [Google Scholar]

[R2] [2].Agračev AA, and Gamkrelidze RV, “ The exponential representation of flows and the chronological calculus,” Mathematics of the USSR-Sbornik, vol. 35, no. 6, pp. 727, 1978. [Google Scholar]

[R3] [3].Albertini F and Sontag ED, “For neural networks, function determines form,” Neural Networks, vol. 6, no. 7, pp. 975–990, 1993. [Google Scholar]

[R4] [4].Bergmann TO, Groppa S, Seeger M, Mölle M, Marshall L, and Siebner HR, “Acute changes in motor cortical excitability during slow oscillatory and constant anodal transcranial direct current stimulation,” Journal of Neurophysiology, vol. 102, no. 4, pp. 2303–2311, 2009. [DOI] [PubMed] [Google Scholar]

[R5] [5].Breakspear M, “ Dynamic models of large-scale brain activity,” Nature Neuroscience, vol. 20, no. 3, pp. 340–352, 2017. [DOI] [PubMed] [Google Scholar]

[R6] [6].Bressan A and Piccoli B, Introduction to the Mathematical Theory of Control, Vol. 1, American Institute of Mathematical Sciences, Springfield, 2007. [Google Scholar]

[R7] [7].Centorrino V, Gokhale A, Davydov A, Russo G, and Bullo F, “Euclidean contractivity of neural networks with symmetric weights,” IEEE Control Systems Letters, vol. 7, pp. 1724–1729, 2023. [Google Scholar]

[R8] [8].Chen R, Singh M, Braver TS, and Ching S, “Dynamical Models Reveal Anatomically Reliable Attractor Landscapes Embedded in Resting State Brain Networks,” Imaging Neuroscience, 2025, doi: 10.1162/imag_a_00442. Available: https://doi.org/10.1162/imag_a_00442. [DOI] [Google Scholar]

[R9] [9].Chen RTQ, Rubanova Y, Bettencourt J, and Duvenaud D, “Neural Ordinary Differential Equations,” Advances in Neural Information Processing Systems. 2018. [Google Scholar]

[R10] [10].Coron JM, Control and nonlinearity. No. 136, American Mathematical Soc, 2007. [Google Scholar]

[R11] [11].Davydov A, Proskurnikov AV and Bullo F (2022, June). “Non-Euclidean contractivity of recurrent neural networks,” In 2022 American Control Conference (ACC) (pp. 1527–1534). IEEE. [Google Scholar]

[R12] [12].DeGiorgio CM, and Krahl SE, “Neurostimulation for drug-resistant epilepsy,” Continuum: Lifelong Learning in Neurology, vol. 19, no. 3, pp. 743–755, 2013. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R13] [13].Forti M and Tesi A, “New conditions for global stability of neural networks with application to linear and quadratic programming problems,” IEEE Transactions on Circuits and Systems I: Fundamental theory and applications, vol. 47, no. 7, pp. 354–366, 1995. [Google Scholar]

[R14] [14].Glorot X, Bordes A, and Bengio Y, “Deep sparse rectifier neural networks,” in Proc. 14th Int. Conf. on Artificial Intelligence and Statistics (AISTATS), pp. 315–323, 2011. [Google Scholar]

[R15] [15].Grover S, Wen W, Viswanathan V, Gill CT, and Reinhart RMG, “Long-lasting, dissociable improvements in working memory and long-term memory in older adults with repetitive neuromodulation,” Nature Neuroscience, vol. 25, no. 9, pp. 1237–1246, 2022. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R16] [16].Gu S, Pasqualetti F, Cieslak M, Telesford QK, Yu AB, Kahn AE, Medaglia JD, Vettel JM, Miller MB, and Grafton ST, and et al. , “Controllability of structural brain networks,” Nature Communications, vol. 6, no. 1, pp. 8414, 2015. [Google Scholar]

[R17] [17].Hirsch MW and Smale S, Differential Equations, Dynamical Systems, and Linear Algebra, Academic Press, New York-London, 1974. Pure and Applied Mathematics, Vol. 60. [Google Scholar]

[R18] [18].Hopfield JJ, “Neurons with graded response have collective computational properties like those of two-state neurons,” Proceedings of the National Academy of Sciences, vol. 81, no. 10, pp. 3088–3092, 1984. [Google Scholar]

[R19] [19].Liu Y-Y, Slotine J-J, and Barabási A-L, “Controllability of complex networks,” Nature, vol. 473, no. 7346, pp. 167–173, 2011. [DOI] [PubMed] [Google Scholar]

[R20] [20].Lohmiller W and Slotine JJE. “On contraction analysis for nonlinear systems,” Automatica, vol. 34 no. 6, pp 683–696, 1998. [Google Scholar]

[R21] [21].Kawski M, Chronological calculus in systems and control theory, Mathematics of Complexity and Dynamical Systems, pp. 88, 2011. [Google Scholar]

[R22] [22].López-Alonso V, Cheeran B, Río-Rodríguez D, and Fernández-del-Olmo M, “Inter-individual variability in response to non-invasive brain stimulation paradigms”, Neuron, vol. 7, no. 3, pp. 372–380, 2014. [Google Scholar]

[R23] [23].Mastrogiuseppe F, and Ostojic S, “Linking connectivity, dynamics, and computations in low-rank recurrent neural networks”, Neuron, vol. 99, no. 3, pp. 609–623, 2018. [DOI] [PubMed] [Google Scholar]

[R24] [24].Medaglia JD, “Clarifying cognitive control and the controllable connectome,” Wiley Interdisciplinary Reviews: Cognitive Science, vol. 10, no. 1, p. e1471, 2019. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R25] [25].Muldoon SF, Pasqualetti F, Gu S, Cieslak M, Grafton ST, Vettel JM, and Bassett DS, “Stimulation-based control of dynamic brain networks”, PLoS Computational Biology, vol. 12, no. 9, pp. e1005076, 2016. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R26] [26].Power JD, Schlaggar BL, Lessov-Schlaggar CN, and Petersen SE, “Evidence for hubs in human functional brain networks”, Neuron, vol. 79, no. 4, pp. 798–813, 2013. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R27] [27].Qin Y, Bassett DS, and Pasqualetti F, “Vibrational control of cluster synchronization: Connections with deep brain stimulation,” in Proc. 61st IEEE Conf. on Decision and Control (CDC), pp. 655–661, 2022. [Google Scholar]

[R28] [28].Qin Y, Nobili AM, Bassett DS, and Pasqualetti F, “Vibrational stabilization of cluster synchronization in oscillator networks,” IEEE Open Journal of Control Systems, vol. 2, pp. 439–453, 2023. [Google Scholar]

[R29] [29].Ramachandran P, Zoph B, and Le QV, “Searching for activation functions,” arXiv preprint, arXiv:1710.05941, 2017. [Google Scholar]

[R30] [30].Sarychev AV, “Lie extensions of nonlinear control systems” Journal of Mathematical Sciences, vol. 135, no. 4, pp. 3195–3223, 2006. [Google Scholar]

[R31] [31].Singh MF, Wang A, Cole M, Ching S, Braver TS, “Enhancing task fMRI preprocessing via individualized model-based filtering of intrinsic activity dynamics”, NeuroImage, vol. 247, pp. 118836, 2022. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R32] [32].Singh MF, Braver TS, Cole MW, and Ching S, “Estimation and validation of individualized dynamic brain models with resting-state fMRI”, NeuroImage, vol. 221, pp. 117046, 2020. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R33] [33].Singh MF, Matthew, Wang A, Anxu, Braver TS, and Ching S, “Scalable surrogate deconvolution for identification of partially-observable systems and brain modeling”, Journal of neural engineering, vol. 17, no. 4, pp. 046025, 2020. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R34] [34].Schiff SJ, Neural control engineering: the emerging intersection between control theory and neuroscience, MIT Press, 2011. [Google Scholar]

[R35] [35].Schuessler F, Dubreuil A, Mastrogiuseppe F, Ostojic S, and Barak O, “Dynamics of random recurrent networks with correlated low-rank structure”, Physical Review Research, vol. 2, no. 1, pp. 013111, 2020. [Google Scholar]

[R36] [36].Sontag ED, and Qiao Y, “Remarks on controllability of recurrent neural networks”, Proceedings of the IEEE Conference on Decision and Control., Vol. 1. Institute of Electrical and Electronics Engineers Inc., 1998. [Google Scholar]

[R37] [37].Tamekue C and Ching S, “Control analysis and synthesis for general control-affine systems,” arXiv preprint, arXiv:2503.05606, 2025. [Google Scholar]

[R38] [38].Tamekue C, Prandi D, and Chitour Y, “A mathematical model of the visual MacKay effect”, SIAM Journal on Applied Dynamical Systems, vol. 23, no. 3, pp. 2138–2178, 2024. [Google Scholar]

[R39] [39].Tucsnak M and Weiss G, Observation and control for operator semigroups, Springer, 2009. [Google Scholar]

[R40] [40].Yi Z, Heng P-A, and Fu AW-C, “Estimate of exponential convergence rate and exponential stability for neural networks”, IEEE Transactions on Neural Networks, vol. 10, no. 6, pp. 1487–1493, 1999. [DOI] [PubMed] [Google Scholar]

[R41] [41].Zhang H, Wang Z, and Liu D, “A comprehensive review of stability analysis of continuous-time recurrent neural networks,” IEEE Transactions on Neural Networks and Learning Systems, vol. 25, no. 7, pp. 1229–1262, 2014. doi: 10.1109/TNNLS.2014.2317880. [DOI] [Google Scholar]

PERMALINK

On the control of recurrent neural networks using constant inputs

Cyprien Tamekue

Ruiqi Chen

ShiNung Ching

Roles

Abstract

I. Introduction

Fig. 1.

II. Representation of solutions

III. Controllability of the neural networks with STEP FUNCTION

A. The case of linear activation functions