We define constrained activation functions φ(z) by properties that mirror those of biological neurons: the boundedness condition is stated as follows:

∃Mφ > 0 s.t. |φ(z)| ≤ Mφ, ∀z ∈ ℝ. (3)

Smoothness and controlled derivatives are given as follows:

φ ∈ C for k ≥ 2, with Lipschitz-continuous derivatives. Specifically,

′ ′′

‖φ ‖∞ ≤ L1, ‖φ ‖∞ ≤ L2, (4)

where ∥⋅∥∞ denotes the supremum norm.

Saturation is as follows:

lim φ (z) = 0. (5)

z→±∞ In contrast, common US activations (e.g., ReLU, φ(z) = max(0, z)) violate (3) and (5) and have undefined or unbounded second derivatives, violating (4).

We define robust neural coding as the stability of the population code fθ(x) under infinitesimal input perturbations.

For an input x and a perturbation vector δ ∈ ℝ , the local sensitivity metric is the amplification factor κ(x, δ): ‖fθ(x + δ) − fθ(x)‖p κ(x, δ) = , ‖δ‖q where ∥⋅∥ and ∥⋅∥q are appropriate p-and q-norms (typically p = q = 2). A network exhibits stable coding in a region ℛ ⊆ i κ(x, δ), where f remains bounded and close to 1 for small ‖δ‖ across all x ∈ ℛ.

The first-order Taylor approximation yields the fundamental link to the network Jacobian:

fθ(x + δ) − fθ(x) ≈ Jf(x)δ, (6)

where,

∂fθ d ×d

J (x) = ∈ ℝ L 0,

Thus, the worst-case local amplification is governed by the operator norm (induced p-norm) of the Jacobian: ∥ Jf(x) δ ∥p κmax(x) = sup =∥ Jf(x) ∥p→q, δ≠0 ∥ δ ∥q For p = q = 2, this is the spectral norm σmax(Jf(x)), the largest singular value of Jf(x).

Therefore, the robustness condition can be formulated as requiring the Jacobian's spectral properties to be controlled:

∥ J(x) ∥2≤ η, ∀x ∈ ℛ, where η ≈ 1 indicates enhanced local stability, and η "1 indicates potential instability.

For a feedforward network, the Jacobian factors into a product of layer-wise Jacobians:

(L) (L) (L−1) (L−1) (1) (1)

Jf(x) = D (x)W D (x)W … D (x)W ,

l ′ (l) ′ (l)

where D (x) = diag (φ (z1 ), … , φ (z )) is a diagonal matrix of activation derivatives at layer l.

dl

This factorization reveals the global sensitivity. ‖Jf(x)‖ is bounded by the product of layer-wise sensitivities:

(l) (l) ′ (l) (l)

‖Jf(x)‖ ≤ ∏ ∥ D (x) ∥⋅∥ W ∥= ∏ (max ∣ φ (z ) ∣) ⋅ σmax(W ).

l=1 l=1

′ (l)

Constrained activations directly enforce maxi ∣ φ (z ) ∣≤ L1, providing a built-in, data-dependent stabilizer that prevents the runaway growth of this product, a mechanism absent in US activations like ReLU, where φ (z) can be identically 1.

We instantiate ensembles of fully-connected networks with depth L ∈ {3,5,8,10} and uniform width dl = n = 128. Parameters are initialized using the initialization scheme [11], scaled to preserve signal variance across layers in expectation: W ∼ (0,2/nin). All networks are analyzed in their untrained, randomly initialized

ij

state to isolate the intrinsic, architecture-driven properties of the nonlinearities.

We systematically compare three families of activation functions φ(z):

The US is stated as follows:

I. Biologically-Constrained, Bounded (BC-B): Hyperbolic Tangent tanh(z), rescaled to have a Lipschitz constant L1 ≈ 1.

II. Biologically-Constrained, Saturating & Smooth (BC-SS): a smoothed, saturated function defined as φβ(z) =

1 βz 1 β(z−1)

log(1 + e ) − log(1 + e ) for β = 4, which approximates a soft threshold with explicit bounds on all derivatives.

All functions are centered and scaled so that z∼ (0,1)[φ(z)] = 0 and Varz∼ (0,1)[φ(z)] = 1 to ensure fair comparison of signal propagation.

For a given input x0 ∼ (0, Id ), we generate a set of perturbations {δi} with norms ϵi =∥

−4, −1

δi ∥2 logarithmically spaced in [10 10 ]. We compute the exact output deviation Δ(ϵi) =∥ fθ(x0 + δi) − fθ(x0) ∥2.

We analyze the scaling law between input and output perturbation magnitudes by fitting the relationship:

log Δ(ϵ) ≈ αlog ϵ + constant, via linear regression over the range of ϵ. The exponent α characterizes robustness:

I. α > 1: superlinear growth (unstable, sensitive).

II. α = 1: linear growth (neutral, Jacobian-dominated).

III. α < 1: sublinear growth (stable, robust).

We report the distribution of α over an ensemble of inputs and network initializations.

We compute the full Jacobian matrix Jf(x) for each input x via efficient backpropagation of the identity matrix. From Jf(x), we extract:

Spectral norm is:

σmax(x) = Jf(x).

r 2 2

Defined as exp(− ∑i=1 pi log pi), where pi = σi / ∑ σj are the normalized squared singular values. This

argument quantifies the dispersion of sensitivity across input directions.

Condition number is:

σmax(x) κ(Jf(x)) = , σmin(x) indicating the anisotropy of the input-output transformation.

We compute the statistics (mean, 95th percentile) of these quantities across the input distribution to assess typical and worst-case sensitivity.

To assess global representational stability, we sample inputs from a simple, low-dimensional latent manifold. Let a latent variable s ∈ [−1,1] parameterize a spiral in input space: x(s) = U(s1, s2, γs1s2, 0, … ,0) + ξ, where U is a random orthogonal matrix, and ξ is small noise. We map 1000 such points through the network to obtain a population activity manifold ℳ = {fθ(x(si))}i=1 .

We quantify manifold stability under input noise ζ ∼ (0, ξ I) using:

I. Distance preservation: the Pearson correlation ρ between pairwise distances in input space, ∥ xi − xj ∥2, and output space, ∥ yi − yj ∥2.

II. Metric distortion: The normalized metric distortion is defined as

1 ∣ ∥ yi − yj ∥2 −∥ xi − xj ∥2 ∣ = ∑ , Np ∥ xi − xj ∥2

where Np is the number of pairs.

III. Dimensionality: the participation ratio Deff = (∑i λi) /(∑ λi ) of the activity covariance matrix, where λi are eigenvalues.

This multi-scale analysis connects microscopic properties (activation derivatives, Jacobian spectra) to macroscopic functional outcomes (perturbation growth, manifold stability). By comparing the UC, BC-B, and BC-SS families under identical architectural and initial conditions, we directly attribute differences in robustness to the imposed mathematical constraints on φ(z).

The perturbation scaling exponent α (defined in Section 3.2) provides a global summary of robustness. For UC ReLU networks, we observe α = 1.02 ± 0.05 (mean ± std across initializations), indicating purely linear, Jacobian-dominated scaling. In stark contrast, networks with bounded, saturating activations exhibit clear sublinear growth (α < 1).

(BC-B - tanh): α = 0.87 ± 0.03.

(BC-SS): α = 0.82 ± 0.04.

This sublinearity implies a compressive nonlinearity at the network level: small input perturbations are damped rather than amplified. The effect intensifies with depth (L). For ReLU networks, α remains near 1 regardless of L. For constrained networks, α decreases monotonically with L (e.g., from ~0.90 at L = 3 to ~0.78 at L = 10 for BC-SS), demonstrating that the stabilizing effect is compositional and accumulates across layers.

The mechanism underlying sublinear growth is directly observable in the spectral properties of the Jacobian Jf(x).

The mean spectral norm x[σmax(Jf] diverges sharply between families. For ReLU networks, σmax scales exponentially with depth in the worst case, with a 95th percentile often exceeding 10 for L = 10 . For constrained networks, σmax remains tightly bound. Crucially, its distribution across inputs is heavily concentrated near 1.

I. UC (ReLU): [σmax] = 4.3 ± 2.1, P95 = 18.7.

II. BC-B (tanh): [σmax] = 1.1 ± 0.2, P95 = 1.8.

III. BC-SS: [σmax] = 0.9 ± 0.1, P95 = 1.3.

The effective rank of Jf(x) reveals how sensitivity is distributed across input directions. ReLU networks exhibit a high effective rank (~70% of the maximum), indicating that perturbation amplification is diffuse across many independent input modes. Constrained networks exhibit a dramatic collapse in sensitivity rank to ~20-30% of the maximum. It indicates that amplification is channeled into a vastly smaller, more structured subset of directions, consistent with the notion of a protected, low-dimensional manifold for population dynamics.

The full singular value spectrum {σi} of Jf reveals the core geometric difference. ReLU networks exhibit a heavy-tailed distribution, with a long "tail" of non-negligible singular values. Constrained networks show a rapid exponential decay of σi. This decay implies that the mapping fθ effectively projects inputs onto a low-dimensional output subspace, thereby drastically attenuating noise components orthogonal to it.

The global consequences of Jacobian stabilization are evident in the structure of population activity manifolds ℳ.

Under input noise (ξ = 0.1), the distance correlation ρ between input and output pairwise distances is significantly higher for constrained networks.

I. UC (ReLU): ρ = 0.65 ± 0.08, = 0.41 ± 0.06.

II. BC-B (tanh): ρ = 0.88 ± 0.04, = 0.18 ± 0.03.

III. BC-SS: ρ = 0.92 ± 0.03, = 0.12 ± 0.02.

The lower distortion for BC-B and BC-SS confirms that the relative geometry of the stimulus set is preserved more faithfully.

The participation ratio Deff of the population activity is consistently 20 − 40% lower for constrained networks compared to their ReLU counterparts of the same width, even when driven by the same 2D latent manifold. It indicates that constrained nonlinearities actively compress representations into lower-dimensional subspaces, thereby inherently regularizing the network's response by discarding extraneous, noise-prone dimensions.

(l) (l)

Deconstructing the Jacobian product Jf = ∏ D W provides the fundamental explanation. The key term is

(l) ′ (l)

the layer-wise gradient matrix. D = diag(φ (z )).

′ ′ (l)

For ReLU, φ (z) ∈ {0,1}. The proportion of active units (φ = 1) varies widely with x, causing ∥ D ∥ to fluctuate between 0 and 1. It can lead to pathological inputs where all layers are fully active, resulting in ∥ Jf ∥≈ ∏ ∥ W ∥,

which grows exponentially with L.

For constrained activations, ∣ φ (z) ∣ is not only bounded above by L1 ≈ 1, but its distribution across units

′ 2 ′ 2

and inputs is strongly sub-Gaussian. We find [φ (z) ] ≪ ( [∣ φ (z) ∣]) Due to saturation for large Jf(x). This means that as signals propagate, the effective Lipschitz constant of each layer, ∥ D ∥It is not just bounded, but typically much less than 1. The product ∏l ∥ D ∥ therefore decays exponentially with depth, directly counteracting the potential growth from weight matrices W . This built-in, data-dependent damping is the origin of the observed sublinear perturbation growth and spectral taming.

A critical concern is that stability might be achieved at the expense of excessive smoothing or loss of feature selectivity. We test this by examining the kernel alignment of each network's infinite-width initial Gaussian Process (GP) kernel K(x, x ) [12] with a simple target kernel (e.g., a Laplace kernel). While ReLU networks induce a kernel with higher absolute variation, constrained networks produce kernels with comparable or superior alignment to smooth target functions. It suggests that the constrained networks are not merely smoothing inputs, but are forming structured, stable representations that can still support complex function approximation. The stability arises from the geometry of the representation, not from a loss of nonlinear expressive power.

Biologically inspired constraints on activation functions' boundedness, saturation, and smoothness serve as a powerful, intrinsic regularizer. They enforce stable sensory coding by: 1) taming the spectral norm of the network Jacobian, 2) inducing sublinear perturbation growth, 3) collapsing the effective sensitivity rank, and 4) promoting low-distortion, geometry-preserving population manifolds. This robustness emerges from the dynamic, data-dependent damping of signal propagation (∥ D ∥< 1), a mechanism inherently absent in US piecewise-linear architectures.