Introduction

Transcranial focused ultrasound (tFUS) has emerged as a promising non-invasive approach for modulating deep brain activity with millimetre-scale spatial precision and access to targets that are difficult to reach with electromagnetic stimulation [5, 6, 17, 27, 38, 49, 50]. Realization of this promise will require addressing the widely observed variability in the direction and magnitude of the brain’s response to stimulation across subjects, sessions, and trials [7, 14, 36, 37, 39, 51, 56]. Potential sources of this variability include biophysical differences between subjects [14], differ-ences in experimental design across studies [37, 56], and the inherently non-stationary nature of neural activity – stimulation interacts with ongoing dynamics, such that the response is itself mod-ulated by the baseline state [39]. In the transcranial magnetic stimulation (TMS) literature, state-dependent effects are well documented [4, 48]. Such findings suggest that explaining response variability will require modeling tFUS as a perturbation of a non-stationary dynamical system.

In practice, most efforts to explain variable tFUS responses have focused on the stimulation itself by mapping dose-response curves across acoustic parameters [36, 56]. Less attention has been given to the choice of outcome measure. The dependent variable is typically a measure of overall activity, such as population firing rate or fluorescence amplitude. These readouts quantify how much activity changed, but are largely agnostic to the underlying temporal organization. In a dynamical system, an amplitude increase can arise from different regimes: a tonic elevation, a transient burst, or a shift in excitatory-inhibitory balance. Moreover, dynamical reorganization and amplitude changes need not co-occur: a perturbation can alter temporal organization without producing a resolvable change in mean amplitude, or vice versa.

This limitation motivates a complementary viewpoint from nonequilibrium statistical physics. The brain is a paradigmatic far-from-equilibrium system that continuously expends energy to main-tain organized activity patterns. tFUS is itself an energetic perturbation: it deposits acoustic energy into dissipative neural tissue. Understanding neuromodulation is thus partly a question of how a system that is already far from thermodynamic equilibrium processes an additional, exogenous energy flux. The field of stochastic thermodynamics [16, 47] provides the framework to address this question. Central to this framework are fluctuation theorems [11, 22] that relate nonequilibrium work to equilibrium free-energy differences, and estimators of entropy production rate (EPR) that quantify broken detailed balance and have revealed nonequilibrium dynamics in living systems [3, 20]. In systems weakly coupled to an environment, stochastic thermodynamics relates fluctua-tions, time-irreversibility, and entropy production, connecting observed trajectories to dissipation and the thermodynamic arrow of time [13, 42, 46]. Recent work has shown that temporal irre-versibility estimated from neural time series captures mesoscopic nonequilibrium structure in the brain, with irreversibility varying across brain states, arousal levels, and pharmacological condi-tions [12, 31, 32, 45]. If sonication changes the nonequilibrium character of neural dynamics (i.e., how irreversibly the system evolves in time) then stochastic thermodynamics may afford a new understanding of ultrasonic neuromodulation.

Here we employ tools from stochastic thermodynamics to conduct a novel analysis of an ex-isting dataset from Murphy et al. [36], which combined fiber photometry of the behaving rodent thalamus with tFUS at varying intensities (Figure 1A,B). We estimated an ordinal surrogate of the entropy production rate (EPR) before, during, and after stimulation (Figure 1C,D) to address the following questions: (i) Does tFUS modulate neural irreversibility either during or immediately after sonication? (ii) Does EPR contain information beyond that of calcium amplitude? (iii) Does prestimulation irreversibility predict the subsequent response to stimulation? We present evidence that all three answers are affirmative and discuss their implications for the mechanistic understand-ing of ultrasonic neuromodulation.

Figure 1.

Results

2.1 Calcium amplitude and irreversibility show distinct dose-response pat-terns

We analyzed fiber photometry recordings from the central medial thalamus (CMT) [36]. Each ana-lyzed trial included a 5 second baseline, 5 seconds of stimulation, and the 10 seconds immediately after tFUS (“recovery”; Figure 1B). Pulsed tFUS (pulse repetition frequency PRF = 2.5 Hz; carrier frequency 550 kHz) was applied either to the recorded region (“on-target”) or to the area 3 mm laterally adjacent to the CMT (“off-target”; Figure 1A). The acoustic intensity spanned six values ranging from 0.3 to 7.4 W/cm². We estimated an ordinal surrogate of the EPR (Figure 1C-D) after pooling trials within each animal-condition (see “Pooled-trial EPR and calcium amplitude mea-sures” in Methods). To evaluate dose-response for both EPR and GCaMP amplitude, we employed mixed-effects models with a common specification: the baseline value was included as a covariate, and the predictors of each model included the condition (on-target versus off-target) and acoustic intensity (both linear and quadratic terms; see Eq. 4).

Calcium amplitude increased monotonically with dose in the on-target condition during both stimulation and recovery, while off-target fluorescence remained low across the intensity range (Figure 2C,D,F). On-target stimulation produced significantly larger calcium responses than off-target in both windows (stimulation: p < 0.001; recovery: p < 0.001). On-target-specific dose structure was significant during stimulation, with both the linear (β_on×I = +7.57, p = 0.006) and quadratic (β_on×I2 = −48.3, p = 0.010) interaction terms reaching significance, consistent with a monotonic increase that saturates at higher intensities. During recovery, calcium amplitude showed only a linear dependence on intensity (β_on×I = +16.34, p < 0.001; quadratic p = 0.746).

Figure 2.

Temporal irreversibility showed a qualitatively different dose dependence. In the on-target condition (Figure 2A,E), both the stimulation and recovery windows showed the largest changes from baseline at moderate intensities (2.3–3.7 W/cm²). Interestingly, the off-target condition (Figure 2B,E) also showed a non-monotonic irreversibility response (albeit milder) despite the absence of a strong fluorescence signal, hinting that sonication can alter dynamics in the absence of an amplitude change.

During stimulation, we did not resolve an on-target main effect for EPR (p = 0.690), nor did we find significant condition-intensity interactions (p = 0.216 and p = 0.872 for the linear and quadratic interactions, respectively). The absence of a significant on-target contrast during stimulation is consistent with the observation that off-target sonication also elevated EPR in this window (Figure 2B,E). In recovery, by contrast, on-target EPR was significantly higher than off-target (p < 0.001), and both on-target-specific dose interaction terms were significant (β_on×I = −0.265, p = 0.002; β_on×I2 = −1.550, p = 0.006). Note that a negative regression coefficient on the quadratic dose term is consistent with an inverted-U response. Overall, the effect of on-target stimulation on EPR emerges immediately after sonication. Table 1 summarises the on-target-specific terms for both readouts and both windows.

Table 1.

EPR dose-response is consistent across animals

The inverted-U EPR dose-response pattern during recovery was clearly visible in 4 of 5 animals, with the peak change from baseline (ΔD_sym) occurring between 1.2 and 3.7 W/cm² (Figure S1). The pattern was not observed with off-target stimulation (Figure S2). Thus, the quadratic group-level effect was not driven by a single strong “responder” and was specific to sonication of the CMT.

No significant GCaMP response at 20 Hz PRF

Murphy et al. [36] reported an offline inhibitory effect when stimulating the CMT with a higher PRF of 20 Hz. We were not able to resolve a main effect of on-target stimulation on calcium amplitude during stimulation (p = 0.676) or recovery (p = 0.821), finding instead that the change in GCaMP remained positive (Figure S3). We discuss potential reasons for this discrepancy in the Discussion. The corresponding analysis of EPR did not resolve a significant main effect of on-target sonication during stimulation (p = 0.940), although a nominally significant EPR reduction was present during recovery (p = 0.044).

2.2 Irreversibility cannot be explained by calcium amplitude alone

Figure 3 depicts the change from baseline for both EPR and GCaMP, shown separately during (panel A) and after (panel B) stimulation. The sharp falloff of irreversiblity at higher acoustic intensities, especially during recovery, suggests that EPR is capturing a distinct aspect of CMT activity. We thus measured the correlation between the two assays of neural response – ΔD_sym and ΔGCaMP – during and immediately after on-target tFUS (Figure 3C-D, respectively).

Figure 3.

During stimulation, the changes from baseline were positively correlated (r = 0.58, p < 0.001), even after conditioning the correlation on intensity (r_partial = 0.63, p < 0.001). During recovery, however, the correlation vanished (r = 0.02, p = 0.927; r_partial = 0.16, p = 0.409), despite the fact that both EPR and GCaMP are significantly different from baseline levels during the recovery period (Table 1).

To test whether recovery EPR (Figure 4B) carries dose information that calcium amplitude cannot account for, we fit a set of mixed models (see “EPR beyond calcium amplitude” in Methods) where the base (reduced) model captures the dependence of EPR on its baseline value and the observed change in GCaMP (Figure 4A). Adding linear and quadratic dose terms significantly increased model fit (on-target stimulation trials only; χ²(2) = 13.55, p = 0.001), indicating the presence of dose structure beyond calcium amplitude. When pooling on-target and off-target trials and adding dose interaction terms, model fit was further improved (χ²(3) = 27.23, p < 0.001). Importantly, the ΔGCaMP covariate was not significant (β = −0.006, p = 0.642) in the full model, indicating that the level of irreversibility during recovery is not associated with the change in calcium amplitude. We calculated a residual EPR by regressing out its baseline and the associated change in GCaMP amplitude. The residualized EPR is shown as a function of dose in Figure 4C. In the on-target condition, the residual is positive from 0.3 to 3.7 W/cm² but negative at 5.4 and 7.4 W/cm². In other words, high-dose sonication produced less irreversibility than expected from its calcium-amplitude response.

Figure 4.

We repeated the hierarchical mixed model analysis for the EPR observed during the stimula-tion period (Figure S4), and again found that dose improved fit beyond baseline EPR the change in GCaMP (on-target only; χ²(2) = 7.28, p = 0.026). As in recovery, the model fit was significantly improved with the addition of on-target versus off-target terms (χ²(3) = 15.27, p = 0.002). How-ever, unlike recovery, the ΔGCaMP covariate was statistically significant (β = +0.181, p < 0.001). In summary, after accounting for baseline value and the change in calcium amplitude, EPR still retains significant dose structure, indicating that sonication reshapes recovery dynamics be-yond what amplitude alone would predict. Complete likelihood-ratio test statistics are provided in Table S1.

2.3 Prestimulation irreversibility predicts stimulation responsiveness

We asked whether the prestimulation level of temporal irreversibility – a single-trial marker of neural state – predicted the magnitude or direction of the response to tFUS. Both EPR and GCaMP responses were computed at the single-trial level for baseline, stimulation, and recovery windows across all six acoustic intensities (n = 408 trials for 2.5 Hz PRF).

Baseline state modulates stimulation EPR

To test whether baseline irreversibility influences the subsequent change in EPR during or after sonication, we partitioned trials into terciles based on their prestimulation EPR. Trials were pooled across intensities and condition such that we con-sidered both on-target and off-target stimulation of the CMT. We then fit a mixed-effects model with baseline-EPR tercile, condition, and their interaction as fixed effects, with a random inter-cept for animal (see “Baseline-state mixed-effects models” in Methods). During stimulation, we found a significant baseline-state × condition interaction (χ²(2) = 11.49, p = 0.003), indicating that the relationship between baseline EPR and the subsequent change differs significantly when stimulating the recorded region (Figure 5A). More specifically, on-target stimulation generated a larger increase in EPR when prestimulation irreversibility was lower. In recovery, the correspond-ing interaction was absent (χ²(2) = 1.05, p = 0.592; Figure 5B), indicating that the baseline-state dependence of EPR is concentrated in the stimulation period rather than persisting into recovery.

Figure 5.

Baseline EPR modulates the change in GCaMP

We also investigated whether prestimula-tion irreversibility predicts the change in GCaMP during and after stimulation (see “Baseline-state mixed-effects models” in Methods). During stimulation, we found a significant interaction be-tween baseline EPR and condition (χ²(2) = 9.19, p = 0.010), indicating that the baseline EPR dependence of the calcium response differed between on-target and off-target tFUS. In recovery, the corresponding interaction was not significant (p = 0.113). Baseline EPR also contributed to overall GCaMP variability beyond condition and intensity, reaching significance in recovery (χ²(2) = 6.45, p = 0.040) and trending similarly during stimulation (χ²(2) = 5.56, p = 0.062).

Figure 6C,D shows the corresponding adjusted means: on-target responses are largest in the low-baseline tercile and reduced in the high-baseline tercile, whereas off-target responses remain near zero across strata.

Figure 6.

To quantify the incremental predictive value of baseline EPR, we fit trial-level on-target-only mixed models in which GCaMP was predicted from baseline GCaMP, acoustic intensity (linear and quadratic), and a random animal intercept (see “Incremental prediction of GCaMP by baseline EPR” in Methods). Adding baseline EPR to that model significantly improved fit during stimulation (χ²(1) = 5.04, p = 0.025; β = −2.78, ΔR² = 0.020) and during recovery (χ²(1) = 4.17, p = 0.041; β = −3.54, ΔR² = 0.017). In both windows, the sign was negative: higher prestimula-tion irreversibility predicted a weaker calcium response after controlling for baseline calcium and dose. Figure 6 shows this relationship via added-variable (partial regression) plots and adjusted tercile means.

Thus, baseline EPR predicts overall response magnitude, with lower-baseline trials showing larger calcium responses on average – a relationship that holds after adjusting for baseline GCaMP and acoustic intensity.

2.4 Sensitivity to ordinal parameters

Because ordinal EPR depends on the embedding dimension m and delay τ, we tested the influence of parameter choices on our findings. For the main analysis presented above, we selected m = 3 to minimize the dimensionality of the forward and reverse probabilities, which increase with m!. The value for the delay τ = 8 was selected from a split-half reliability sweep of baseline irreversibility across τ = {2, 4, 6, 8}.

The core findings were preserved with (m, τ) = (3, 6). During recovery, EPR dose-response showed a significantly negative quadratic dose relationship (β_I2 = −1.90, p = 0.032). The EPR during recovery was shown to contain information beyond that of calcium amplitude (on-target-specific dose: χ²(3) = 17.50, p = 5.6 × 10⁻⁴; on-target-only dose: χ²(2) = 8.75, p = 0.013; ΔGCaMP covariate remained null: p = 0.992). Baseline EPR significantly influenced the subse-quent change in both EPR (baseline-EPR × condition interaction: χ²(2) = 7.00, p = 0.030) and GCaMP amplitude (χ²(2) = 6.70, p = 0.035) during stimulation.

The main claims were also largely preserved with a higher-dimensional embedding, namely (m, τ) = (4, 8). The EPR during recovery exhibited a significant quadratic dose dependence (β_I2 = −2.12, p = 0.015), and both on-target-specific recovery dose terms remained significant (β_on×I = −0.39, p = 0.002; β_on×I2 = −2.12, p = 0.015). Moreover, calcium amplitude could not explain the EPR measured during recovery (on-target-specific dose: χ²(3) = 24.59, p = 1.9 × 10⁻⁵; on-target-only dose: χ²(2) = 11.07, p = 0.0040; ΔGCaMP covariate again remained null: p = 0.911). Baseline EPR predicted the subsequent change in irreversibility during sonica-tion in a target-specific manner (χ²(2) = 12.59, p = 0.0018). Baseline irreversibility tended to also predict the associated change in calcium amplitude during sonication, although the effect did not reach significance (p = 0.105).

Discussion

By employing tools from stochastic thermodynamics, this study found that low-intensity focused ultrasound perturbs brain activity in a manner that is not captured by conventional amplitude-based readouts. The GCaMP signal recorded via fiber photometry reflects aggregate calcium dynamics across the illuminated neural population, encompassing both somatic and neuropil compartments [9, 26]. When stimulating the CMT with tFUS, the amplitude of the GCaMP signal increased monotonically with the intensity of the acoustic stimulus [36] (Figure 2). By contrast, the EPR showed a non-monotonic dose-response that peaked at moderate intensity (Figure 2A,B,E), a pat-tern that was visible in both on-target and off-target conditions.

An inverted-U dose–response for temporal irreversibility

The inverted-U pattern of EPR is reminiscent of nonlinear dose-response relationships that have been documented across several forms of brain stimulation and, more broadly, across biology. The classical Yerkes–Dodson law posits that behavioural performance peaks at intermediate levels of arousal or stimulus intensity [52], and the hormesis literature documents biphasic dose-response curves in which low doses stimulate and high doses inhibit across a wide range of biological measures [8]. In transcranial direct-current stimulation, moderate current densities produce the strongest effects on cortical ex-citability and cognition, whereas higher intensities can reverse or abolish these effects [2, 15]. Deep brain stimulation of the central thalamus likewise enhances working memory at low currents and impairs it at higher currents [33]. The concept of stochastic resonance – where an intermedi-ate level of noise optimally enhances signal detection in a nonlinear system – provides a formal framework for understanding such inverted-U phenomena [19, 35]. The findings here point to the possibility that the brain’s response to exogenous energy input is fundamentally nonlinear, with an intermediate regime in which stimulation most effectively engages endogenous dynamics.

For ultrasound, one interpretation is that moderate acoustic doses perturb the system enough to displace it from its resting attractor without overwhelming the endogenous dynamics that gen-erate temporal irreversibility. At higher intensities, the exogenous drive may dominate the local dynamics – producing large calcium transients – while constraining the system to a more stereo-typed, less dissipative trajectory. In the language of nonequilibrium physics, moderate sonication may optimally couple to the system’s internal degrees of freedom, generating the greatest depar-ture from time-reversal symmetry. By contrast, strong sonication may effectively “entrain” the population [21, 41], increasing the amplitude of the response while reducing its dynamical com-plexity. Although speculative, this framing makes the testable prediction that the peak of the EPR dose-response should shift with parameters that alter the balance between exogenous drive and endogenous dynamics, such as pulse repetition frequency or duty cycle.

The post-sonication dissociation between irreversibility and amplitude

It is notable that the dissociation between GCaMP and EPR is most pronounced immediately after sonication (Fig-ure 3B,D; Figure 4). During this period, the system is reorganising after perturbation, and both measures remain elevated above their respective baselines. The fact that the two measures are uncorrelated during this window (Figure 3D) indicates that amplitude and irreversibility reflect distinct aspects of the region’s relaxation. In other words, the path by which the system relaxes after perturbation carries information that is invisible to the amplitude readout. Moderate-intensity sonication may push the system into a region of state space from which the return to baseline in-volves the greatest thermodynamic cost – the largest entropy production — whereas high-intensity stimulation, despite producing larger fluorescence changes, may push the system along a more direct, less dissipative relaxation path.

Off-target stimulation elevates irreversibility without a calcium response

An interesting as-pect of the present findings is that stimulating a region 3 mm adjacent to the site of the recording fiber produced a clear inverted-U EPR elevation during stimulation that was visible in multiple animals (Figures 2B, S2). In the mouse brain, a 3 mm offset from the CMT places the acoustic focus well outside of the thalamus proper [40], yet the two sites remain within reach of common thalam-ocortical and corticothalamic projection systems. It is therefore plausible that off-target sonication perturbs polysynaptically connected circuits whose altered dynamics propagate to the recorded CMT population. If these network-level perturbations produce changes in temporal structure with-out substantially increasing the aggregate calcium flux, they would appear in the irreversibility measure but not in the bulk fluorescence signal. This interpretation is consistent with the broader notion that EPR is sensitive to changes in the temporal organization of activity that need not man-ifest as amplitude changes. Note that for conventional neuromodulation studies, off-target stim-ulation serves as an experimental control. Here, it represents an informative finding in of itself, suggesting that sonication can alter the nonequilibrium character of neural dynamics at a distance, and that EPR may be sensitive to low-amplitude perturbations.

Baseline irreversibility as a susceptibility marker

Another aspect of our findings suggests that EPR serves as a proxy for neural “susceptibility”: trials with lower baseline irreversibility showed systematically higher responses during stimulation, measured with both GCaMP and EPR itself (Figures 5, 6). This is conceptually aligned with findings in other stimulation modalities showing that prestimulus brain state governs the response to transcranial stimulation. In TMS, for exam-ple, cortical excitability tracks the instantaneous phase of endogenous oscillations [44], and state-dependent protocols can selectively induce plasticity that is absent when stimulation is delivered irrespective of brain state [4, 58]. The EPR effect – lower baseline D_sym predicting larger ΔD_sym – could in principle reflect regression to the mean, since a high-baseline trial will mechanically tend toward a negative change score. However, our data argue against such an artifactual explana-tion. The relationship between baseline irreversibility and the subsequent change was significantly modulated by whether tFUS was applied on-or off-target (Figure 5A). Moreover, baseline EPR predicted the GCaMP amplitude response even after controlling for baseline GCaMP and acoustic intensity (Figure 6). This cross-modal prediction (i.e., prestimulation irreversibility anticipates the magnitude of the calcium-amplitude response) is evidence that baseline EPR captures a genuine axis of stimulation susceptibility.

No effects when stimulating at 20 Hz

Murphy et al. [36] reported that 20 Hz PRF stimulation of the CMT produced an offline, inhibitory calcium response. While our analysis of 2.5 Hz stim-ulation was highly congruent with that of Murphy et al. [36], we did not find significant changes in either calcium amplitude or EPR at a PRF of 20 Hz (Supplementary Figure S3). This discrepancy may stem from differences in the preprocessing that were employed in the respective studies, par-ticularly the approaches used to correct for motion artifacts. The finding that 20 Hz stimulation did not modulate irreversibility is consistent with the interpretation that the 20 Hz protocol elicited an overall weaker neuromodulatory effect compared to 2.5 Hz.

Limitations: thermodynamic interpretation

Our employment of an ordinal surrogate of EPR (D_sym; see Eq. 2) draws on the formal connection between time-reversal asymmetry and dissi-pation in nonequilibrium statistical physics. The foundational result – that the Kullback–Leibler divergence between forward and reverse trajectory probabilities equals the mean entropy produc-tion – was established for Markovian systems coupled to a thermal bath [11, 22, 47] and extended to discrete-state systems [16, 42, 43]. Practical estimators that detect broken detailed balance from time-series data have extended the applicability of this framework [34, 46]. Broken detailed bal-ance has been identified at mesoscopic scales in active biological matter [3, 20]. Applications to neural data have shown that human brain dynamics violate detailed balance at large spatial scales, with the degree of irreversibility tracking arousal state and cognitive load [12, 31, 32, 45].

However, neural tissue violates the assumptions of these theorems in several ways. The GCaMP signal reflects calcium dynamics filtered through an indicator with slow kinetics. The system has multiple chemical reservoirs, and the ordinal embedding captures only a low-dimensional projec-tion of the underlying state space. The symmetrised KL divergence used here is therefore best understood as a nonparametric measure of time-arrow asymmetry in the observed signal: a quan-tity that is related to, but not identical with, the thermodynamic EPR of the underlying biophysical system. Nevertheless, our results indicate that this time-asymmetry measure carries information about the neural response to sonication that is inaccessible to conventional amplitude readouts.

Limitations: EPR estimation

EPR estimation from noisy biological time series is intrinsically challenging. Any KL-based irreversibility estimator depends on the embedding parameters and on finite-sample bias, especially in short windows. We employed an ordinal-based method because it is nonparametric, insensitive to monotone rescalings of the signal, and naturally suited to short segments [1, 53]. Alternative approaches to estimating entropy production from time series include visibility-graph irreversibility [24], neural-network classifiers trained on the arrow of time [46], model-based estimators that infer dissipation from partially observed Langevin dynamics [34], and current-based methods that quantify dissipation from fluctuating probability currents [28]; a recent review provides a systematic comparison [54]. Analytical treatments of EPR in neural-field models have also begun to connect irreversibility to network parameters [30].

Broader implications

Nonequilibrium statistical physics offers a language for asking how an exogenous stimulus does work on a living system that is already far from thermodynamic equilib-rium. In this framework, sonication does not merely increase or decrease neural activity; it can also reshape the path which the system takes during and after perturbation. The present results suggest that moderate doses may drive the system into a regime with greater temporal irreversibility than would be expected from fluorescence amplitude alone, whereas higher doses can produce large calcium increases with less irreversibility. Future work could test this more directly by combin-ing time-asymmetry markers with electrophysiology and metabolic imaging to establish whether irreversibility changes co-occur with shifts in biophysical dissipation, and by testing closed-loop stimulation designs that target different baseline states [29, 57]. More generally, our findings sug-gest that stochastic thermodynamics serves as a valuable complement to conventional amplitude-based readouts in neuromodulation, as it captures the interaction between exogenous energy and the brain’s ongoing nonequilibrium dynamics.

Methods

4.1 Dataset

We analysed fiber photometry recordings from the publicly available dataset of Murphy et al. [36], which investigated the neuromodulatory effects of transcranial focused ultrasound (tFUS) in the mouse deep brain. Data was accessed via https://dataverse.harvard.edu/dataset. xhtml?persistentId=doi:10.7910/DVN/PCWRAD. In brief, adult C57BL/6J mice were injected with adeno-associated viral vectors encoding the genetically encoded calcium indicator GCaMP6s [9] under the CaMKII promoter and implanted with 200 µm-diameter optical fibers targeting the thalamic central medial nucleus (CMT). Ultrasound was delivered through the intact skull using a wearable, miniaturised ring transducer operating at a fundamental frequency of 550 kHz with a lateral focal spot of ∼2.3 mm (full-width at half-maximum). The transducer was chronically mounted on a 3D-printed head frame, enabling repeated stimulation sessions in freely moving animals.

Dual-wavelength fiber photometry was performed using interleaved 470 nm (calcium-dependent GCaMP6s excitation) and 405 nm (calcium-independent isosbestic) illumination, providing a noise reference for correcting motion and hemodynamic artifacts. Data was sampled at 32 Hz. Each recording session comprised multiple trials at six output levels corresponding to reported spatial-peak pulse-average acoustic intensities of 0.3, 1.2, 2.3, 3.7, 5.4, and 7.4 W/cm². A trial consisted of a 5-second prestimulation baseline, a stimulation epoch (5 s when stimulating at 2.5 Hz PRF; 40 s at 20 Hz PRF), and a post-stimulation recovery period such that the total length of each trial was 180 s. For our analysis, we retained only the first 10 second of the recovery period.

Stimulation was conducted on-target (sonication directed at the recorded CMT region) and separately off-target (sonication of the area 3 mm laterally adjacent to the CMT; Figure 1A). For each configuration, two pulse repetition frequencies (PRFs) were tested: 2.5 Hz and 20 Hz. Murphy et al. [36] reported a strong dose-dependent GCaMP increase with 2.5 Hz and a milder offline inhibition at 20 Hz. The data analyzed in the present study comprised 5 on-target and 5 off-target animals at 2.5 Hz (210 on-target and 198 off-target trials total) and 5 on-target and 4 off-target animals at 20 Hz (216 on-target and 168 off-target trials total). Note that we excluded one animal from analysis (“cmtb3”) as it did not show any GCaMP response during or after stimulation (Supplementary Figures S6–S7).

4.2 Preprocessing

Each trial was preprocessed independently using a three-step procedure. First, the isosbestic (405 nm) reference channel was regressed out of the GCaMP signal using ordinary least-squares regression fit to the prestimulation baseline period only, such that we did not inadvertently remove stimulation-evoked responses. The regression included an intercept term, and the fitted model was applied to the full trial to remove motion and hemodynamic artifacts common to both channels. Next, a linear trend was fit to the full trial-length residual and subtracted to remove slow drift. Finally, the detrended signal was z-scored using the mean and standard deviation of the prestim-ulation baseline period, yielding a baseline-normalized signal in units of standard deviations. No temporal smoothing was applied, preserving the full temporal structure of the signal for subsequent irreversibility analysis. Note that our preprocessing differs somewhat from that of Murphy et al. [36], which is relevant when comparing findings across studies.

4.3 Ordinal entropy production rate estimation

To quantify temporal irreversibility in the preprocessed calcium signal, we employed ordinal pat-tern analysis [1]. This approach maps a scalar time series into a sequence of discrete symbols by replacing each delay-embedded vector with its ordinal pattern (the rank ordering of its compo-nents). The distribution of ordinal patterns carries information about the temporal structure of the signal, and its asymmetry under time reversal provides a nonparametric estimator of the entropy production rate (EPR) [42, 53, 55].

Specifically, for a time series segment {x_t} of length N, we constructed delay embeddings of dimension m and delay τ:

Each embedding vector v_t was mapped to one of m! possible ordinal patterns by recording the rank ordering of its components. The same procedure was applied to the time-reversed signal {x_N+1−t}, yielding a backward ordinal pattern distribution. The symmetrized Kullback–Leibler (KL) divergence between the forward and backward distributions served as our EPR estimator:

where P_fwd and P_bwd denote the normalized ordinal pattern frequency distributions of the forward and reversed signal, respectively, and where the KL divergence between two discrete distributions P and Q over the m! ordinal patterns π_i is:

Additive Laplace smoothing (ε = 10⁻⁶ per bin) was applied before normalization to avoid undefined logarithms.

Symmetrized versus strict EPR estimator

The trajectory-level EPR is rigorously defined as the one-sided KL divergence between forward and time-reversed path distributions, D_fwd ≡ D_KL(P_fwd∥P_bwd) [42, 47]. We adopted the symmetrized form D_sym (the Jeffreys divergence [23]) following the ordinal-pattern irreversibility literature [53, 55], which employs it due to its improved finite-sample stability: it spreads the cost of any zero-frequency pattern symmetrically across both directions, rather than placing it asymmetrically on whichever distribution happens to lack support at a given element. To verify that this choice does not influence our reported inferences, we re-ran the dose-response analysis of Figure 2 using D_fwd in place of D_sym. The resulting irreversibility measures D_sym and D_fwd were essentially proportional (r > 0.99; median D_sym/D_fwd ratio ≈ 2.0). This is expected from the small-divergence regime in which our ordinal-pattern esti-mates lie: to leading order in ∥P_fwd − P_bwd∥, both KL directions reduce to the same functional, so D_KL(P_fwd∥P_bwd) ≈ D_KL(P_bwd∥P_fwd) and D_sym ≈ 2 D_fwd [10, Ch. 11]. We therefore retained D_sym as our irreversibility measure throughout.

Parameter selection

We employed m = 3 and τ = 8 samples as the parameters for measuring ordinal irreversibility. We fixed m = 3 to keep the ordinal alphabet at 6 states, which promotes stable estimates compared to higher-dimensional embeddings. The delay parameter τ was then se-lected outcome-blind from the set {2, 4, 6, 8} using a baseline-only reliability sweep at fixed m = 3. The delay value τ = 8 exhibited the highest split-half correlation for animal-level baseline EPR and was therefore adopted as the primary setting. Sensitivity analyses were conducted at neighbouring delay scales and with a higher-dimensional embedding to assess robustness to alternative choices (see Section 2.4).

Temporal windowing

Each trial was segmented into three temporal windows: a prestimulation baseline (5 s, the maximum prestimulation time recorded by Murphy et al. [36]), stimulation (5 seconds for 2.5 Hz PRF; 40 s for 20 Hz PRF), and post-stimulation recovery (10 s beginning immediately after sonication ends). A 0.25-second boundary trim was applied at the baseline–stimulation transition to avoid edge effects.

4.4 Statistical analysis

We employed both animal-and trial-level analyses depending on the question being addressed. For the dose-reponse and EPR versus GCaMP analyses of Figures 2, 3, and 4, irreversibility was estimated after pooling trials of each animal and forming one estimate per animal. To investigate the influence of baseline irreversibility on the subsequent response to stimulation (Figures 5,6), EPR was estimated separately for each trial to afford greater statistical power, trading off sample size with the SNR of each individual estimate.

Pooled-trial EPR and calcium amplitude measures

For each (animal × intensity × condition × window) cell, we computed one D_sym value by pooling ordinal-pattern counts across all trials belonging to that cell and then evaluating Eq. 2 on the pooled distribution (see Section 4.3). For GCaMP, the dependent variable was the mean fluorescence amplitude averaged across trials within the same cell. The changes from baseline were defined as ΔD_sym = D_sym^window − D_sym^baseline and ΔGCaMP = GCaMP^window −GCaMP^baseline for EPR and calcium amplitude, respectively. We note that change scores were computed primarily for visualization of effects – when conducting formal inference (see next section), we opted for ANCOVA-style mixed models where baseline value was included as a covariate.

Mixed-effects models

For formal inference, we fit linear mixed-effects models (LMMs) with animal as a random intercept [25], using the same ANCOVA-style specification for both EPR and GCaMP:

where j indexes the animal, i indexes the repeated measure (one sample per intensity × condition), y denotes either D_sym or mean GCaMP, the superscript “window” denotes the temporal window (stimulation or recovery), I_c is the centered acoustic intensity, OnTarget_{i j} is a binary indicator equal to 1 for on-target and 0 for off-target, and u _j ∼ N (0, σ ²) is the random intercept for an-imal j. Including the baseline value as a covariate adjusts for variability in prestimulation state while avoiding the mechanical baseline–change coupling that arises from using difference scores as the outcome. The main-effect dose terms (β_I, β_I2) capture any shared dose dependence across conditions, while the interaction terms β_on×I and β_on×I2 test whether the linear and quadratic dose dependence is specific to on-target stimulation. Models were fit using restricted maximum likeli-hood (REML) with the L-BFGS solver as implemented in the statsmodels Python package.

EPR beyond calcium amplitude

To test whether irreversibility carries dose information beyond what calcium amplitude can explain, we fit a hierarchy of three nested mixed-effects models on the pooled-trial data. All three shared a random intercept for animal.

Model A (baseline + amplitude only) included only prestimulation EPR and the window-matched fluorescence change as predictors:

where ΔGCaMP_{z,i j} is the standardized change in mean calcium amplitude for repeated measure i of animal j. Model B (shared dose) added linear and quadratic intensity terms that were common to both on-and off-target conditions:

Model C (condition-specific dose) added an on-target fixed effect and further allowed the dose dependence to differ between on-target and off-target stimulation:

We employed likelihood-ratio (LR) tests to assess significant increases in explained variance: comparing Model B to Model A tests whether acoustic dose explains any EPR variance beyond baseline state and calcium amplitude change. Comparing Model C to Model B tests whether the dose–EPR relationship differs between on-target and off-target stimulation. Finally, to ask whether dose shapes EPR within the on-target condition alone, we refit Models A and B on the on-target subset only and compared them with a third LR test.

4.5 Baseline neural state and response prediction

To test whether prestimulation temporal irreversibility predicts the response to sonication, we mea-sured baseline EPR at the level of individual trials, focusing the analysis on the 2.5 Hz PRF record-ings.

Trial-level EPR and GCaMP computation

For each trial, we computed ordinal EPR and mean GCaMP amplitude separately for the baseline, stimulation, and recovery windows. Trial-level change scores were defined as ΔD_sym = D_sym^window − D_sym^baseline and ΔGCaMP = ^window − ^baseline, where denotes the mean calcium amplitude within the specified window.

Baseline-state mixed-effects models

To test whether baseline state modulates the response to tFUS, we pooled trials across acoustic intensities and discretised baseline EPR into terciles (Low, Mid, High). For the EPR outcome measure, we fit a hierarchy of mixed-effects models according to:

where α_{Ii j} denotes fixed effects of acoustic intensity (entered as a categorical factor with six levels), S_{i j} ∈ {Low, Mid, High} is the baseline-EPR tercile, γ_{Si j} is the corresponding tercile fixed effect, δ_{Si j} captures the tercile-by-condition interaction, and u _j ∼ N (0, σ ²) is a random intercept for animal j. LR tests compared the Main model to the Null (baseline-state main effect) and the Interaction model to the Main (baseline-state × condition interaction).

The choice of a categorical form for intensity reflects the non-monotonic dose-response of EPR (Figure 4B) – modeling intensity as categorical makes no assumption about the shape of the dose-response and thus avoids introducing misspecification into the estimate of the quantity-of-interest, namely the baseline-on-target interaction.

For GCaMP, the same three-model hierarchy was applied with calcium amplitude as the out-come and baseline GCaMP added as a covariate to adjust for prestimulation fluorescence:

Incremental prediction of GCaMP by baseline EPR

To quantify whether baseline EPR car-ries predictive information about the subsequent GCaMP response beyond what baseline GCaMP and acoustic dose already explain, we fit trial-level mixed models on the on-target trials only. The reduced model predicted GCaMP from its baseline value, centred intensity (linear and quadratic), and a random animal intercept:

The full model included baseline D_sym as an additional predictor:

A significant LR test comparing the full model to the reduced model indicates that baseline EPR explains GCaMP variance that is not captured by baseline fluorescence or dose. Both models were fit with ML estimation. The added-variable (partial regression) plots in Figure 6A,B visualise the residual association between baseline EPR and post-window GCaMP after removing the shared covariates [18].

Supplementary Material

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Additional information

Funding

HHS | National Institutes of Health (NIH) (5R16GM145496)

• Jacek P Dmochowski