PINN-ing the Balloon: A Physically Informed Neural Network Modelling the Nonlinear Haemodynamic Response Function in Functional Magnetic Resonance Imaging

Accurate characterisation of the haemodynamic response function (HRF) is central to interpreting blood-oxygen-level-dependent (BOLD) signals in functional magnetic resonance imaging, yet standard estimation approaches remain centred around phenomenological formulations lacking biophysical grounding. We present a physics-informed neural network (PINN) framework that bridges these paradigms by embedding the Balloon-Windkessel model directly into the training objective of a multi-headed Neural Network. Our aproach simultaneously estimates probable latent neurovascular state variables such as cerebral blood inflow, metabolic rate of oxygen consumption, blood volume, and deoxyhaemoglobin content, through an indirect optimisation scheme in which the predicted BOLD signal is obtained via convolution of the estimated HRF with experimental stimuli. Training is governed by a composite loss, balancing differential-equation residuals, physiological initial conditions and data fidelity. In simulations with temporal signal-to-noise ratios representative of clinical acquisitions, the framework recovered ground-truth state variables with coefficients of determination exceeding 0.99 and mean squared errors below 10 −3, at a physics-to-data weighting of 0.40:0.60. Application to 1.5 T block-design fMRI data from an ischaemic stroke patient yielded physiologically plausible, subject-specific HRF estimates, establishing feasibility of single-subject, physics-constrained HRF inference without reliance on fixed gamma basis assumptions.To our knowledge, this constitutes the first deployment of a single PINN incorporating the full Balloon-Windkessel model within an indirect training objective, reconstructing full BOLD observations, positioning PINN-based haemodynamic modelling as a principled and personalised route towards more interpretable and patient-specific fMRI biomarkers. 1/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint

Physically Informed Neural Networks; Haemodynamic Response Function; fMRI; Balloon model; Differential equations. 1 Introduction 1 Magnetic Resonance Imaging (MRI) is a cornerstone of clinical and research 2 neuroimaging; as it captures both structural and physiological data, besides being 3 non-invasive. Functional MRI (fMRI) enables the study of brain function in vivo by 4 measuring fluctuations in the blood oxygenation level-dependent (BOLD) signal, which 5 arise from haemodynamic and metabolic responses to neuronal activity whether task or 6 state related [43]. 7 The haemodynamic response function (HRF; Fig. 1) is an idealised, noise-free 8 representation of the BOLD response to a brief stimulus (typically < 4 s) [39,40,44]. 9 Commonly, the HRF is thought to reflect the brain physiological response, a 10 macroscopic signature of neurovascular coupling shaped by vascular and metabolic 11 dynamics [12]. A canonical HRF model arises from studies of haemodynamic responses 12 in healthy individuals and is typically estimated as a combination of two gamma13 probability density functions, where the first models the shape of the initial 14 stimulus-response, while the second models the undershoot [20,44]. 15 HRF models are crucial for identifying brain activation, and it is often desirable to 16 estimate their parameters with physiological interpretations [50], particularly when 17 exploring the relationship between experimental stimuli and brain responses [10]. 18 However, the neural mechanisms linked to metabolism and blood vessel dynamics that 19 generate the BOLD signal remain an active area of research [5,38]. 20 Various methods and techniques have been proposed to estimate the HRF e.g. 21 Poisson-based approaches [21], Gaussian parameterisations [48], and combinations of 22 gamma functions [52]; Bayesian non-parametric estimation [13] and BOLD 23 deconvolution methods that exploit stimulus timing [26,29,50]. 24 Several toolboxes support BOLD simulation, e.g., The Virtual Brain, Dynamic 25 Causal Modelling (DCM) in SPM, and NeuRosim [38]. In brief, The Virtual Brain 26 typically generates BOLD-like signals by convolving a neural mass/mean-field activity 27 time course with a canonical HRF, whereas DCM uses variants of the Balloon model, a 28 widely used generative model of cerebral haemodynamics [6,7,42,54].29 The Balloon model, initially proposed by Buxton et al. [7] (see Section 2.1), 30 represents the venous compartment of a region of interest (ROI) as an expandable 31 “balloon”. The resulting nonlinear dynamical system is commonly written in terms of 32 four normalised state variables that capture key physiological processes: cerebral blood 33 inflowf in(t), cerebral metabolic rate of oxygen consumption m(t), cerebral blood 34 volume v(t), and deoxyhaemoglobin content q(t) [21,50]. A primary challenge, however, 35 lies in precisely identifying this dynamic system, primarily due to noise in experimental 36 data [50]. 37 During the last years, a new type of artificial neural networks that leverages the 38 purely data driven methods by using physical equations that describes the process 39 under study was developed. This Physics-informed neural networks (PINNs), integrate 40 observational data with mechanistic knowledge by embedding governing equations 41 directly into the training objective [34,41]. The embedded equations constrain learning 42 to physically admissible solutions, support the solution of forward and inverse problems, 43 and can improve robustness when data are sparse or noisy [47]. By leveraging automatic 44 differentiation, PINNs compute derivatives without finite-difference discretisation,45 which can be advantageous for parameter inference in dynamical models.46 2/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint From a Machine Learning (ML) perspective, incorporating prior knowledge is a 47 powerful strategy for tackling key challenges, including limited training data, improving 48 model generalisation, and ensuring the physical plausibility of results [47]. From a 49 numerical methods perspective, solving differential equations using artificial neural 50 networks (ANNs) differs significantly from classical numerical methods. The solution to 51 the equation comes from optimising the network parameters, which do not depend on 52 the equation’s dimensions, thereby mitigating the problems and drawbacks associated 53 with dimensionality [60]. 54 In this work, we develop a PINN-based approach for HRF estimation by combining 55 BOLD observations with the haemodynamic constraints of the Balloon model [6,7,21]. 56 The method leverages the flexibility of neural networks while enforcing biophysical 57 structure on blood flow, volume, and oxygenation dynamics. In this way, it aims to 58 bridge purely data-driven HRF fitting and fully specified generative haemodynamic 59 modelling, yielding BOLD reconstructions that remain physiologically interpretable. 60 Figure 1. HRF descriptors: (1) Peak Amplitude (Height HP), (2) Time to Peak (TTP), (3) Full Width at Half Maximum (FWHM), (4) Time to Onset (TO), (5) and area under the curve of the first peak (AUC). If there is an undershoot: (6) Minimum Undershoot Height (MU), (7) Time to Undershoot Minimum, (8) and time to recover baseline value (TT0). 2 Material and methods 61 2.1 The Balloon Model 62 The Balloon model expresses the BOLD signalh(t) (Eq. 1) as a static, nonlinear 63 function of blood volume v(t) and deoxyhaemoglobin content q(t). 64 h(t) = V0 ·  k1 · [1 − q(t)] + k2 ·  1 − q(t) v(t)  +k 3 ·[1− v(t)]  k1 = 4.3 V0 ϑ0 E0 TE k2 = ϵ r0 E0 TE k3 = 1 − ϵ Here V0 is the resting blood volume fraction; ϑ0 is the frequency offset at the outer 65 surface of a magnetised vessel for fully deoxygenated blood at 1.5 T; E0 is the resting 66 oxygen extraction fraction; TE is the echo time; ϵ is the ratio of intra- and extravascular 67 3/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint Figure 2. Balloon model with parallel driven CBF and CM RO2, coupling the increase in q with the normalised CM RO2 (m) directly (see [6]), which is driven by a second input signal ( ICM RO2) like the normalised CBF (increased to 0.05). Further, fout is described by the equations of [6], which causes v to decrease more slowly. signals; andr 0 is the slope relating the intravascular relaxation rateR ∗ 2I to oxygen68 saturation. Although alternative parameter values and model extensions can be found 69 in [6,15,42,54], we adopt the values of k1, k2, and k3 reviewed by Stephan et al. [54], 70 which have been used in subsequent work [11,14,24,30,31,38,45]. 71 As shown in Eq. 1, and given a mean transit time τM T T, the dynamics of v(t) are 72 driven by inflow fin(t) through an outflow term fout(v, t). The deoxyhaemoglobin state 73 q(t) is jointly shaped by blood volume, outflow, and the metabolic drive m(t) [7]. 74 dq dt = 1 τM T T ·[m(t)− q(t) v(t) ·f out(v, t)] dv dt = 1 τM T T ·[f in(t)−f out(v, t)] fout(v) = v 1 α (t) +τ dv dt (1) Here α is the Grubb stiffness exponent (a proxy for venous compliance) [28,54], and 75 τcontrols the time scale of the viscoelastic outflow adjustment.76 We follow the two-input pipeline in Fig. 2 proposed by [38], which builds on 77 extensions suggested by [21] and [5]. The key modelling choice is to separate the inflow 78 fin(t) from the metabolic drive m(t), allowing them to be driven in parallel by the same 79 stimulus input I(t) through Eq. 2. 80 d2fin dt2 +κ· d f dt =λ f ·I(t) − γ · (fin(t) − 1) d2m dt2 +κ· dm dt =λ m ·I(t) − γ · (m(t) − 1) (2) We choose λf such that fin(t) exhibits an underdamped response (allowing overshoot 81 and undershoot), whereas λm is chosen to yield a critically damped response for m(t). 82 The parameters κ and γ control signal decay and feedback regulation, respectively. 83 2.2 Physically Informed Neural Network 84 A standard PINN architecture comprises an input layer, one or more hidden layers, and 85 an output layer. The input is a feature vector sampled from the problem domain; the86 hidden layers learn a latent representation; and the output layer produces the final 87 prediction vector [37,46]. 88 Within this framework, we represent the unknown solutions of Eqs. 2 and 1 with a 89 single deep neural network uθ(t), where θ denotes all trainable parameters (weights and 90 4/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint biases). We implement a multi-headed (parallel-output) multilayer perceptron, adapted 91 from [46], to simultaneously produce the four Balloon state variables. 92 Our model (Fig. 3) takes a uniformly sampled time vector as input (30 s, sampled 93 every centisecond), uses two hidden layers (128 → 256 units, 256 → 512 units), and 94 applies SoftMax activations in the shared trunk. The output layer comprises four linear 95 “heads” (Softplus activation with β = 1.0 and threshold = 20.0), corresponding to the 96 Balloon state variables: ˆ uf and ˆ um, estimates of fin(t) and m(t) from Eq. 2, and ˆ uv 97 and ˆ uq, estimates of v(t) and q(t) from Eq. 1. The coupling among heads mirrors the 98 model structure where fin(t) drives v(t), and (v (t), m(t)) jointly determine q(t). The 99 network output is a 4 × 3000 array. 100 The HRF is not a direct network output. Instead, it is obtained by substituting ˆ uv 101 and ˆ uq into Eq. 1, and the resulting h(t) is then used to reconstruct the BOLD signal 102 via Eq. 7. 103 Training used random initiail value for θ according to Glorot initialisation [25] and 104 ADAM optimiser with 10000 iterations, an initial learning rate of 1 × 10−3, and a 105 step-wise decay factor of 0.15 every 1000 iterations. 106 All experiments were implemented in Python 3.10.12 using PyTorch 2.3.1 and CUDA 107 12.1. Computations were performed in single precision floating-point representation on 108 a single NVIDIA GeForce RTX 4080 GPU, on 13th Gen Intel Core i9-13900 x 32. 109 Ubuntu 24.04.4 LTS, GNU-Linux 6.8.0-106-generic as OS. 110 The manuscript was written using the free version of PrismAI, an AI-powered online 111 LaTeX editor. Grammarly was used as a paid web browser add-on during manuscript 112 writing. The free version of Claude AI Sonnet 4.6 was used for proofreading paragraphs 113 and optimising homemade Python code during the analysis of experiments. 114 Figure 3. PINN Architecture: We used a three-layer network architecture. The input layer processes the main input from a time array through a softMax activation function. Then a two-stage hidden layer, also connected through a softmax, ends up feeding into the output layer. The output layer uses a softPlus activation function, obtaining m(t) and f(t). In addition to the hidden layers, v(t) and q(t) are also fed with the outputs of the fin and m(t) heads. Blues indicate the layers, redish lines indicate the activation function. Orange circles indicate the output functions A central advantage of PINNs is their independence from the dimensionality of the 115 equations to be solved, as training can be formulated directly as an optimisation 116 problem i.e, the identification of a set of the network’s parameters ( ˆθ) that minimises a 117 total loss function Ltot (eq 3). When training for differential equations, Ltot is defined 118 as a weighted sum of the equation residual term Leq, the initial-condition term Lic, and 119 the data-fitting term Ldata, following Eq. 3. 120 5/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint Ltot =ω eqLeq +ω icLic +ω dataLdata (3) For the theoretical term Leq, representing the unknown solution of the Balloon 121 model by a deep neural network, allows us to define the ODEs residuals as 122 Rθ(˚t) = duθ d˚t d˚t dt − G[u θ] (˚t) whereGis a differential operator, ˚tthe normalised time domain; thus, minimising 123 the residuals leads to a theoretical (physical) component of the loss function 124 parametrised by ˚t. For further details on the standardisation of variables and the 125 reparameterisation of equations, see 5. Then Leq can be expressed as 126 Leq = 1 Nt NtX i≥I0 |Rθ(˚ti)|2 (4) In our case, the Balloon model is formulated as a system of ordinary differential 127 equations; thus, {ti}Nt i=1 denotes the set of standardised time-domain samples. 128 To solve Eq. 2, we apply a single brief impulse (a 1 s boxcar). Its onset time tI0 129 marks the end of the resting baseline period used to impose initial conditions. 130 Concretely, we constrain not only the states at t = 0 (i.e., fin(0), m(0), v(0), q(0)) but 131 also all samples with t < t I0 (i.e., fin(t < t I0), m(t < t I0), v(t < t I0), q(t < t I0)), which 132 encourages the network to maintain physiologically plausible basal levels before 133 stimulation. While the forward Balloon model is commonly solved with Dirichlet initial 134 conditions, we additionally impose Cauchy-type conditions (state and derivative 135 constraints) to reduce spurious oscillations around the baseline. 136 In Eq. 5, the first term penalises deviations of the network output from the baseline 137 state, and the second term penalises deviations of its first derivative from the baseline 138 derivative, both evaluated for t < t I0. The quantities x0 and x0 denote the baseline 139 values for the state variables and their derivatives, respectively. 140 Lic = 1 Ni NiX i<I0 |uθ(ti)−x 0|2 +| duθ dt (ti)−x 0|2 (5) To define the data misfit term Ldata, we use an indirect training objective. The 141 predicted BOLD signal ˆY (t) (Eq. 6) is obtained by convolving the estimated HRF ( ˆh(˚t)) 142 with the stimulus function s(t) (Eq. 7), and is fitted to either simulated or in vivo data. 143 The HRF ˆh(˚t) is computed by substituting the network estimates ˆ uv(˚t) and ˆ uq(˚t) for 144 v(t) and q(t) into Eq. 1. 145 Ldata = 1 Nd NdX i=1 | ˆY(t i)−y (ti)|2 (6) ˆY (t) = s(t) ∗ ˆh(˚t) (7) Here {y(ti)}Nd i=1 denotes the observed data samples indexed by time. 146 Therefore, the PINN must solve the dynamical system in Eqs. 2 and 1 while 147 simultaneously explaining the observed BOLD signal through the forward observation 148 model in Eqs. 1 and 7. 149 6/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint 2.3 Experiments 150 Our aim is to assess PINN’s ability to estimate the HRF by regressing the convolution151 of the balloon model’s output with known stimuli against an observed BOLD signal. 152 This requires, before the convolution, solving the coupled ODE systems in Eqs. 1 and 2 153 We evaluate the performance under three experimental settings. First, retrieving the 154 HRF from noiseless BOLD Simulated Data. Secondly, after we can make sure our155 proposal is able to restor the Balloon model state variables under perfect conditions, we 156 commited to retrieve the HRF from noisy BOLD Simulated Data: In these experiments, 157 we trained our PINN, repeating the previous paradigm, introducing additive Gaussian 158 noise before subsampling to emulate due tSNR and TR condition; after that we face 159 naturally the application of our framework to real BOLD data from a block-designed 160 fMRI study previously performed: We proceed to train the PINN against in vivo data. 161 Each experiment was repeated 100 times; runs yielding implausible HRF descriptors 162 were discarded, considering empirical criteria similar to those defined by [52] and used 163 by [9]:0e.g peaks occurring during or before the impulse function I from Eqs. 2, BOLD 164 peaks signal of magnitud greater then 15%, T T P > 15 s. 165 2.3.1 Simulated data 166 In the first two experiments, we train the PINN using simulated data generated through 167 home-programmed of the Balloon model. We simulate the response to a single 1 s 168 impulse I(t) in Eq. 2, propagate the resulting states through Eq. 1, to compute h(t) 169 using Eq. 1, and finally generate BOLD time series by convolving h(t) with a 288 s 170 stimulus train, equivalent to an experimental protocol comprising ten 3 s blocks 171 separated by irreguar intervals. Parameter values for Eqs. 1, 2 and 1 follow [6,38] and 172 are summarised in Table 4. 173 For the second simulation setting, following [52], we add independent Gaussian noise 174 to achieve temporal SNR (tSN R ) equivalent to in-vivo conditions 175 Y(t) = Y (t) + (b+ σe) e ∼ N(0, 1) (8) Here Y(t) corresponds to the noise-added simulated fMRI data;Y (t) denotes the 176 ideal noiseless BOLD simulation from the convolution, while b = 116.8 corresponds to 177 the base level and σ = 0.0135 to the additive noise standard deviation in order to 178 emulate the real tSN R ∼ 70 of our in vivo data. After simulation and noise addition, 179 we subsampled Y(t) to mimic 165 volumes (T R = 1.75 s; approximately 5 min 180 acquisition). 181 2.4 Application to Real Data 182 We next apply the model to observational data to assess practical HRF recovery and 183 interpretability. Data come from a study on acute ischaemic stroke conducted between 184 March 2022 and September 2023, approved by the Regional Ethics Committee under 185 Resolution N o.200 − 2026, and carried out following the 2013 revision of the Declaration 186 of Helsinki. The patient, a 52-year-old male, presented with an ischaemic core in the 187 right thalamic region measuring 0.55 ml in volume, with a National Institutes of Health 188 Stroke Scale (NIHSS) score of 1 and an evolution time of 7 hours before hospital arrival. 189 Images were obtained using an eight-element head coil on a 1.5T Signa HDxt 190 scanner (General Electric, Milwaukee, WI, USA). For fMRI, T2∗-weighted echo-planar 191 imaging (EPI) was employed with a repetition time (TR) of 1.75 s, echo time (TE) of 192 60 ms, and a spatial resolution of 1 .9 × 1.9 × 5mm3. The stimulation paradigm involved 193 7/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint passive bilateral wrist flexion-extension movements [2], manually performed by a third 194 party at approximately 1 Hz. Movement timing was coordinated using visual cues on a 195 screen, signalling the start and stop of each motion. Each wrist movement lasted 3 196 seconds, with inter-stimulus intervals ranging from 15 to 30 seconds, yielding 11 197 event-related activations. The total acquisition time for the fMRI session was 5 minutes. 198 After standard preprocessing, slice timing correction was applied, followed by 199 realignment (estimation and reslicing), co-registration of functional and anatomical 200 images, and smoothing with a full width at half maximum (FWHM) of 6 mm in each 201 direction (6, 6, 6) using SPM12 software (Wellcome Trust Centre for Neuroimaging, 202 London, UK). The postcentral gyri were defined as the regions of interest (ROI) for 203 BOLD signal extraction. The ischaemic lesion side of the patient determined the 204 designation of the ischaemic postcentral gyrus (PCG) and the non-ischaemic PCG. We 205 applied baseline correction by subtracting the mean BOLD signal value from the whole 206 time series. 207 Due to random weight initialisation of the PINN, each experiment was performed 208 100 times, yielding 100 HRF estimates ˆh(t). Following [9], we characterise each ˆh(t) 209 using the descriptors shown in Fig. 1: height to peak (HP), time to peak (TTP), 210 FWHM, time to onset (TO; time to a 10% increase from baseline), and the area under 211 the first peak (AUC; using trapezoidal rule). When an undershoot is present, we also 212 report the minimum undershoot height (MUH), time to undershoot minimum (TUM), 213 and time to return to baseline (TT0). After visual inspection, we validated by 214 identifying implausible descriptor values: any value that fell outside predefined 215 acceptable ranges was considered invalid, its estimated HRF was deemed abnormal, and 216 it was therefore excluded from subsequent analysis; after which, mean and confidence 217 intervals (CI) for the estimations and their descriptors were calculated using 218 bootstrapping with 10000 iterations. Implausible values included but are not limited to: 219 nought HP, TTP or FWHM greater than 10 s, TO lower than 2 s. 220 3 Results 221 3.1 Noiseless Simulation 222 The PINN model required approximately 2000 consecutive training iterations to 223 converge within a stable total loos function, but 10000 iterations were used in each 224 training. After ruling out ∼ 10% of ˆh, the time series estimated by PINN for each state 225 variable are illustrated in Fig. 5, within its standard deviation. We compare the 226 normalised blood inflow, cerebral metabolic rate of oxygen, volume, and 227 deoxyhaemoglobin concentration time series predicted by the PINN model against the 228 numerical solutions (our “ground truth” ), finding reasonable agreement between them 229 when measured using Coefficient of determination ( R2), mean squared error (MSE), L2 230 relative error (L 2RE)and Spearman Correlation (Sρ) (Table 1). We observe that our 231 network closely approximates the numerical state variables, while a comparison of its 232 HRF descriptors (shown in Fig. 1), is shown in Fig. 4 233 For our experiment on BOLD noiseless simulated data, we obtained an MSE lower 234 than 10−3. The L2RE between our estimations and the ground truth signal was 235 ∼ 10−3. The MSE ∼ 10−5 between the reconstruction and the data. These elements, 236 along with other goodness-of-fit metrics, are summarised in the table 1 237 3.2 Noise added Simulation 238 Again, the total loss function required approximately 1800 consecutive training 239 iterations to get stable but 10000 iterations were used in each training. Having trained 240 8/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint Figure 4. HRF estimation descriptors (Noiseless simulation). In blue, the descriptions for the ground truth, in green, the description of a double gamma estimation. The box plots show the diversity of solutions given by the PINN: in orange, the median; the red dot shows the mean; and the red line shows the standard deviation. Figure 5. Balloon model state variable and BOLD reconstruction from noiseless simulated data, with standar deviation of the valid training outputs. From left to right: fin and m from Eq. 2, Balloon Model core: v and q from Eq. 1, HRF from Eq. 1, BOLD estimation from the convolution of the hrf estimation and boxcar stimulus of 3 s as shown un 7. Second row: BOLD signal reconstruction by convolution with 3-s stimuli. In all the plots, dashed lines show the numerical solution, whereas continuous lines show the output from ⃗uθ, vertical lines show the impulse/stimulus onset. R2 RMSE (10−3) L2RE(10 −3) Sρ f(t) 0.9882 5.066 0.4979 0.6881 m(t) 0.9546 2.488 2.478 0.6380 v(t) 0.9554 1.126 1.120 0.9984 q(t) 0.9995 0.353 0.355 0.9998 hrf(t) 0.9998 0.029 13.11 0.9855 Table 1. Goodness of fit metrics for a mean regression against noiseless simulation. R2 the coefficient of determination; RMSE root mean squared error; L2RE relative error and Sρr Spearman Correlation coefficient 9/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint Figure 6. HRF estimation descriptors for Noise added simulation ( tSN R = 70.1). In blue, the descriptions for the ground truth. The box plots show the diversity of solutions produced by the PINN’s several training runs: in orange, the median; the red dot shows the mean of our results. R2 RMSE (10−3) L2RE(10 −3) Sρ f pinn 0.996 2.833 2.785 0.948 m pinn 0.995 0.843 0.839 0.877 v pinn 0.994 0.426 0.426 0.999 q pinn 0.995 1.181 1.187 0.980 hrf predict 0.996 0.130 0.130 0.957 Table 2. Goodness of fit metrics for a regression against simulation with white noise added, tSNR = 70.1. R2 the coefficient of determination; RMSE root mean squared error; L2RE relative error and Sρr Spearman Correlation coefficient the PINN model, and after ruling out ∼ 10% of ˆh, using the same criteria as in the 241 previous experiment, the time series estimated by PINN for each state variable is 242 illustrated in Fig. 7, next to standard deviation. As before, we compare the estimated 243 state variables against the numerical “ground truth”, finding reasonable agreement 244 between them (Table 2). We observe that, with weight ratio of 04:06 between physics 245 and data, our PINN model is capable of closely approximating the numerical state 246 variables even after noise was added. 247 3.3 Application to Real Data 248 As seen in the left panel of Fig. 8 and 9, for our experiment estimating the HRF from 249 fMRI data, we got a stable overall MSE around 10 0 after approximately 2000 training 250 iteration; a decomposition of the total loss in its weighted addends can be seen in the 251 right panel of both figures. We allowed the iterations on each training run to go past 252 4000, as it was seen that some networks’ initial set of parameters could converge 253 belatedly 254 Inter hemisferic comparison and description of the estimated HRF can me seen in 255 Fig. 11. Here, except for TO all other descriptors have non-overlapping interquartile 256 ranges, sugesting the we are actually capturing two distinct elements: more prolonged 257 haemodynamic response (FWHM), a larger AUC and deeper undershoot consistent with 258 a higher peak response magnitude (HP) finding on the right side. The right hemisphere 259 takes longer to reach its undershoot minimum (TTU) and to return to baseline(TT0), 260 suggesting slower haemodynamic recovery overall. 261 10/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint Figure 7. Balloon model state variable and BOLD reconstruction from simulated data with noise addition (tSNR=70.1), with the standard deviation from the valid trainings. From left to right: fin and m from Eq. 2, core of the Balloon Model ( v and q) from Eq. 1, HRF from Eq. 1, BOLD estimation given a boxcar stimulus of 3s as shown in 7. In all the plots, dashed lines show the numerical solution, whereas continuous lines show the output from ⃗uθ. The vertical lines mark the beginning of the impulse/stimulation. Figure 8. Loss function example of a training run against the ischaemic right hemisphere BOLD data: in the left panel, the total loss, weighted sum of the elements on the right panel; the loss from the data, the loss from Cauchy’s initial conditions, the loss from the Dirichlet initial conditions, and the loss from the differential equations. Along with the BOLD reconstruction, A comparison between the hemispheric latent 262 state variables can be seen in Fig. 10 sugesting that these differeces could be due to 263 variations in the deaoxyhaemogline dynamic (q (t)) 264 Given that we have no access to a ground truth (state variables) for these in vivo 265 data; Table 3 shows fitness metrics regarding the BOLD signal reconstruction 266 moderated fiting and showing positive correlation . 267 4 Discussion 268 Estimation of the HRF from fMRI is a non-invasive method to investigate the brain’s 269 haemodynamic response to task-driven neuronal activation. Phenomenological 270 approaches [13,22,26,40,52] have been helpful so far; however, they lack a suitable 271 biophysical correlate. 272 11/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint Figure 9. Loss function example of a training run on the non-ischaemic left hemisphere BOLD data: in the left panel, the total loss, weighted sum of the elements on the right panel; the loss from the data, the loss from Cauchy’s initial conditions, the loss from the Dirichlet initial conditions, and the loss from the differential equations. R2 RMSE L2RE Sρ Left BOLD(t) 0.405 1.008 0.771 0.663 Right BOLD(t) 0.503 1.068 0.705 0.770 Table 3. Goodness of fit metrics for a regression against in vivo data, left(right) non-ischemic (ischemic) hemisphere BOLD signal, with tSNR=66.2 (76.3). R2 the coefficient of determination; RMSE root mean squared error; L2RE relative error and SρrSpearman Correlation coefficient Physics-informed neural networks (PINNs), sometimes also referred to as 273 theory-trained neural networks, incorporate mechanistic constraints (e.g., differential 274 equations) directly into neural network training [37]. Introduced by [47], PINNs leverage 275 the universal approximation capabilities of ANNs while penalising violations of the 276 assumed governing equations, making them well suited for solving differential equations 277 (forward problem) and performing parameter inference (backward problem) [32,34,37].278 Our PINN approach to HRF estimation using block design fMRI data, can be cast 279 as an optimisation problem that does not require prior knowledge of gamma probability 280 distributions or suitable priors to Bayesian inference; instead, we used the Balloon 281 model equations, constraining the regression to avoid overfitting to noisy data, while 282 also providing a plausible ad hoc biophysical explanation from which to estimate 283 haemodynamic responses using in vivo clinical data. 284 This approach relies on a balance between the weights assigned to physical 285 knowledge and data. In our work, we show that for simulated signals with tSNR values 286 similar to those of in vivo data acquired by our team, with a standard 1.5 T scanner 287 under clinical conditions, we are able to retrieve the theoretical ground truth with a 288 proportion (wode/wdata) of 0 .40 : 0.60. When the same proportion is applied to clinical 289 data, the retrieved state variables deviate from the vanilla solution (forward problem) of 290 the balloon model equations, hinting at elements of a dynamic specific to the patient 291 from whom data were collected. Real fMRI data contains physiological noise, motion 292 artefacts, and neural variability that simulated data does not, our results show that the 293 PINN is not overfitting to noise The right hemisphere has a higher tSNR (76.3) than 294 the left (66.2). Better signal quality on the ischaemic side provides the PINN with a 295 cleaner target, which likely contributes to the superior goodness of fit metrics observed 296 there.The present work was a proof-of-concept effort. 297 The clinical applications of the present paradigm are not immediate; first of all, the 298 12/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint Figure 10. Balloon model state variable and BOLD reconstruction for the ischemic right hemisphere BOLD signal (blue, tSNR=76.3),the left non-ischemic hemisphere BOLD signal (red, tSNR=66.2) with the standard deviation from valid training runs. First row, from left to right: fin and m from Eq. 2, core of the Balloon Model ( v and q) from Eq. 1, HRF from Eq. 1, Second and third rows, BOLD estimation given boxcar stimuli of 3s as shown in 7. In all the plots, dashed lines show the numerical solution, whereas continuous lines show the output from ⃗uθ. The vertical lines mark the beginning of the impulse/stimulation. Figure 11. HRF estimation descriptors comparison for the ischemic right hemisphere BOLD signal (blue, tSNR=76.3), the left non-ischemic hemisphere BOLD signal (red, tSNR=66.2). The box plots show the diversity of solutions produced by the PINN’s several training runs. The star marker signals the descriptive value for the theoreti- cal/numeric reference, and the dot marker shows the mean of our results. The band around the HRF estimation shows the standard deviation from valid training runs 13/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint training data came from an specific experimental design, and even when PINNs can be 299 used to perform parameter estimation, inverse problems are known for multiple 300 solutions; thus our result constitute one of many sets of state variables reconstructing 301 the BOLD signal given our parametrisation, regularisation or variations of loss functions, 302 its components or scalar weights determination might help with this endeavour. 303 The results we present here were obtained under computational conditions that 304 support iterative training and repeated model evaluation in a short time, resources that 305 may not be routinely available in standard clinical environments. 306 5 Conclusions 307 When constrained to solve the Balloon Windkessel equations, our PINN-based 308 framework, in combination with fMRI data at a conventional magnetic field (1.5 T) and 309 TR(1.75 s), estimates not only the haemodynamic response function but also a plausible 310 origin based on the model’s state variables. This estimation method produces 311 physiologically consistent and personalised results, as it can be trained on single-subject 312 data. To our knowledge, this is the first study to apply a single PINN for a highly313 nonlinear differential equation system with the inclusion of the bold equation and the 314 convolution in an indirect training. Beyond model refinement, future work could 315 broaden its applicability by including a more extensive participant pool and diverse 316 fMRI paradigms, while improving robustness through sensitivity analyses and 317 modifications to the network architecture. 318 Author Contributions 319 Rodrigo H. Avaria: Conceptualisation, Data Curation, Formal Analysis, Investigation, 320 Methodology, Software, Validation, Visualization, Writing − Original Draft Preparation; 321 David Ortiz: Methodology, Writing − Review & Editing; Javier Palma−Espinosa: 322 Methodology, Visualization, Writing − reviewing and editing; Astrid Cancino: Data 323 Curation, Writing − Review & Editing; Pablo Cox: Data Curation, Writing − Review 324 & Editing; Rodrigo Salas: Conceptualisation, Funding Acquisition, Supervision, Writing 325 − Review & Editing; Steren Chabert: Conceptualisation, Funding Acquisition, Project 326 Administration, Resources, Supervision, Writing − Review & Editing. 327 Acknowledgments 328 The authors would like to acknowledge funding from the National Agency for Research 329 and Development (ANID) through FONDECYT N ◦1231268, and N ◦3250439, and 330 ANID-Millennium Science Initiative Programme ICN2021 004.331 Financial Disclosure 332 None reported. 333 Conflicts of Interest 334 The authors declare no conflicts of interest. 335 14/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint

336 1. C. Altamura, M. Reinhard, M.-S. Vry, C. P. Kaller, F. Hamzei, F. Vernieri, P. M. 337 Rossini, A. Hetzel, C. Weiller, and D. Saur. The longitudinal changes of bold 338 response and cerebral haemodynamics from acute to subacute stroke. a fmri and 339 tcd study.BMC Neuroscience, 10(1):151, Dec. 2009. 340 2. T. Arichi, G. Fagiolo, M. Varela, A. Melendez-Calderon, A. Allievi, N. Merchant, 341 N. Tusor, S. J. Counsell, E. Burdet, C. F. Beckmann, and A. D. Edwards. 342 Development of bold signal haemodynamic responses in the human brain. 343 NeuroImage, 63(2):663–673, Nov. 2012. 344 3. J. T. Ashburner, S. J. Kiebel, T. E. Nichols, and W. D. Penny, editors. Statistical 345 Parametric Mapping. Elsevier Science, Burlington, 2011. Description based upon 346 print version of record. 347 4. B. Bonakdarpour, T. Parrish, and C. Thompson. haemodynamic response 348 function in patients with stroke-induced aphasia: Implications for fmri data 349 analysis. NeuroImage, 36(2):322–331, June 2007. 350 5. R. B. Buxton. The thermodynamics of thinking: connections between neural 351 activity, energy metabolism and blood flow. Philosophical Transactions of the 352 Royal Society B: Biological Sciences, 376(1815):20190624, Nov. 2020. 353 6. R. B. Buxton, K. Uludaˇ g, D. J. Dubowitz, and T. T. Liu. Modeling the 354 haemodynamic response to brain activation. NeuroImage, 23, 2004. 355 7. R. B. Buxton, E. C. Wong, and L. R. Frank. Dynamics of blood flow and 356 oxygenation changes during brain activation: the balloon model. Magnetic 357 resonance in medicine, 39(6):855–864, 1998. 358 8. V. Calhoun, T. Adal?, M. Kraut, and G. Pearlson. A weighted least-squares 359 algorithm for estimation and visualization of relative latencies in event-related 360 functional mri. Magnetic Resonance in Medicine, 44(6):947–954, 2000. 361 9. S. Chabert, R. Salas, E. Cantor, A. Veloz, A. Cancino, M. Gonz´ alez, F. Torres, 362 and C. Bennett. haemodynamic response function description in patients with 363 glioma. Journal of Neuroradiology, 51(4):101156, June 2024. 364 10. W. Chen. haemodynamic Response Function Modeling. PhD thesis, University of 365 North Carolina at ChapelHill, 2012. 366 11. Y. Chen, Q. Bukhari, T. W. Lin, and T. J. Sejnowski. Functional connectivity of 367 fMRI using differential covariance predicts structural connectivity and behavioral 368 reaction times.Network Neuroscience, 6(2):614–633, June 2022. 369 12. F. Cignetti, E. Salvia, J.-L. Anton, M.-H. Grosbras, and C. Assaiante. Pros and 370 cons of using the informed basis set to account for haemodynamic response 371 variability with developmental data. Frontiers in Neuroscience, 10, July 2016. 372 13. P. Ciuciu, J.-B. Poline, G. Marrelec, J. Idier, C. Pallier, and H. Benali. 373 Unsupervised robust nonparametric estimation of the haemodynamic response 374 function for any fmri experiment. IEEE Transactions on Medical Imaging, 375 22(10):1235–1251, Oct. 2003. 376 14. C. Coronel-Oliveros, C. Gießing, V. Medel, R. Cofr´ e, and P. Orio. Whole-brain 377 modeling explains the context-dependent effects of cholinergic neuromodulation. 378 NeuroImage, 265:23, Jan. 2023. 379 15/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint 15. T. L. Davis, K. K. Kwong, R. M. Weisskoff, and B. R. Rosen. Calibrated 380 functional mri: Mapping the dynamics of oxidative metabolism 381 (hypercapniacerebrovascular reactivity). Neurobiology, 95:1834–1839, 2 1998. 382 16. S. Deguchi and M. Asai. Dynamic & norm-based weights to normalize imbalance 383 in back-propagated gradients of physics-informed neural networks. Journal of 384 Physics Communications, 7(7):19, July 2023. 385 17. J.-R. Duann, T.-P. Jung, W.-J. Kuo, T.-C. Yeh, S. Makeig, J.-C. Hsieh, and T. J. 386 Sejnowski. Single-trial variability in event-related bold signals.NeuroImage,387 15(4):823–835, Apr. 2002. 388 18. J. V. Duarte, C. Guerra, C. Moreno, L. Gomes, and M. Castelo-Branco. Changes 389 in haemodynamic response function components reveal specific changes in 390 neurovascular coupling in type 2 diabetes. Frontiers in Physiology, 13, Jan. 2023. 391 19. A. Frau-Pascual, F. Forbes, and P. Ciuciu. Physiological models comparison for 392 the analysis of asl fmri data. Proceedings - International Symposium on 393 Biomedical Imaging, 2015-July:1348–1351, 7 2015. 394 20. K. Friston, P. Fletcher, O. Josephs, A. Holmes, M. Rugg, and R. Turner. 395 Event-Related fMRI: Characterizing Differential Responses. NeuroImage, 396 7(1):30–40, Jan. 1998. 397 21. K. Friston, A. Mechelli, R. Turner, and C. Price. Nonlinear responses in fmri: 398 The balloon model, volterra kernels, and other haemodynamics. NeuroImage, 399 12:466–477, 10 2000. 400 22. K. J. Friston, editor. Statistical parametric mapping: the analysis of funtional 401 brain images. Elsevier/Academic Press, Amsterdam ; Boston, 1st ed edition, 2007. 402 23. K. J. Friston, N. Trujillo-Barreto, and J. Daunizeau. Dem: A variational 403 treatment of dynamic systems. NeuroImage, 41:849–885, 7 2008. 404 24. M. Gatica, F. E. Rosas, P. A. M. Mediano, I. Diez, S. P. Swinnen, P. Orio, 405 R. Cofr´ e, and J. M. Cortes. High-order functional redundancy in ageing 406 explained via alterations in the connectome in a whole-brain model. PLoS 407 Computational Biology, 18(9):21, Sept. 2022. 408 25. X. Glorot and Y. Bengio. Understanding the difficulty of training deep 409 feedforward neural networks. In Y. W. Teh and M. Titterington, editors, 410 Proceedings of the Thirteenth International Conference on Artificial Intelligence 411 and Statistics, volume 9 of Proceedings of Machine Learning Research, pages 412 249–256, Chia Laguna Resort, Sardinia, Italy, 13–15 May 2010. PMLR. 413 26. G. H. Glover. Deconvolution of impulse response in event-related bold fmri1. 414 NeuroImage, 9(4):416–429, Apr. 1999. 415 27. C. Grady. The cognitive neuroscience of ageing. Nature Reviews Neuroscience, 416 13(7):491–505, June 2012. 417 28. R. L. GRUBB, M. E. RAICHLE, J. O. EICHLING, and M. M. 418 TER-POGOSSIAN. The Effects of Changes in PaCO2 Cerebral Blood Volume, 419 Blood Flow, and Vascular Mean Transit Time. Stroke, 5(5):630–639, Sept. 1974. 420 Publisher: American Heart Association. 421 16/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint 29. L. Hernandez-Garcia and M. O. Ulfarsson. Neuronal event detection in fmri time 422 series using iterative deconvolution techniques. Magnetic Resonance Imaging, 423 29(3):353–364, Apr. 2011. 424 30. R. Herzog, P. A. M. Mediano, F. E. Rosas, P. Lodder, R. Carhart-Harris, Y. S. 425 Perl, E. Tagliazucchi, and R. Cofre. A whole-brain model of the neural entropy 426 increase elicited by psychedelic drugs. Scientific Reports, 13:12, Apr. 2023. 427 31. R. Herzog, P. A. M. Mediano, F. E. Rosas, A. I. Luppi, Y. S. Perl, 428 E. Tagliazucchi, M. Kringelbach, R. Cofr´ e, and G. Deco. Neural mass modelling 429 for the masses: Democratising access to whole-brain biophysical modelling with 430 FastDMF, Apr. 2022. 431 32. S. Huang, W. Feng, C. Tang, Z. He, C. Yu, and J. Lv. Partial Differential 432 Equations Meet Deep Neural Networks: A Survey. IEEE Transactions on Neural 433 Networks and Learning Systems, 36(8):13649–13669, Aug. 2025. 434 33. J. Jimenez and T. Ozaki. Local linearization filters for non linear 435 continuous-discrete state space models with multiplicative noise. International 436 Journal of Control, 76:1159–1170, 8 2010. 437 34. G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang. 438 Physics-informed machine learning. Nature Reviews Physics, 3(6):422–440, May 439 2021. 440 35. I. E. Lagaris, A. Likas, and D. I. Fotiadis. Artificial neural networks for solving 441 ordinary and partial differential equations. IEEE Transactions on Neural 442 Networks, 9:987–1000, 1998. 443 36. S.-K. Lee, H. W. Yoon, J.-Y. Chung, M.-S. Song, and H. Park. Analysis of 444 functional mri data based on an estimation of the haemodynamic response in the 445 human brain. Journal of Neuroscience Methods , 139(1):91–98, Oct. 2004. 446 37. K. Luo, J. Zhao, Y. Wang, J. Li, J. Wen, J. Liang, H. Soekmadji, and S. Liao. 447 Physics-informed neural networks for PDE problems: a comprehensive review. 448 Artificial Intelligence Review, 58(10):323, July 2025. 449 38. O. Maith, H. ¨Ulo Dinkelbach, J. Baladron, J. Vitay, and F. H. Hamker. Bold 450 monitoring in the neural simulator annarchy. Frontiers in Neuroinformatics, 16, 3 451 2022. 452 39. G. Marrelec, H. Benali, P. Ciuciu, and J.-B. Poline. Bayesian estimation of the 453 haemodynamic response function in functional mri. In AIP Conference 454 Proceedings, volume 617, pages 229–247. AIP, 2002. 455 40. G. Marrelec, H. Benali, P. Ciuciu, M. P´ el´ egrini-Issac, and J. Poline. Robust 456 bayesian estimation of the haemodynamic response function in event-related bold 457 fmri using basic physiological information. Human Brain Mapping, 19(1):1–17, 458 Mar. 2003. 459 41. C. Meng, S. Seo, D. Cao, S. Griesemer, and Y. Liu. When physics meets machine 460 learning: A survey of physics-informed machine learning. arXiv preprint 461 arXiv:2203.16797, Mar. 2022. 462 42. T. Obata, T. T. Liu, K. L. Miller, W. M. Luh, E. C. Wong, L. R. Frank, and 463 R. B. Buxton. Discrepancies between bold and flow dynamics in primary and 464 supplementary motor areas: application of the balloon model to the 465 interpretation of bold transients. NeuroImage, 21:144–153, 1 2004. 466 17/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint 43. S. Ogawa and Y.-W. Sung. Functional magnetic resonance imaging. Scholarpedia, 467 2(10):3105, 2007. 468 44. F. Pedregosa-Izquierdo. Feature extraction and supervised learning on fMRI: from 469 practice to theory. PhD thesis, Universit´ e Pierre et Marie Curie - Paris VI, 2016. 470 45. F. Protopapa, S. Kulashekhar, M. J. Hayashi, R. Kanai, and D. Bueti. Effective 471 connectivity in a duration selective cortico-cerebellar network. Scientific Reports, 472 13(1):17, Nov. 2023. Publisher: Nature Publishing Group. 473 46. P. Protopapas and H. Raval. Generalized one-shot transfer learning of linear 474 ordinary and partial differential equations. In The Symbiosis of Deep Learning 475 and Differential Equations III, page 7, 2023. 476 47. M. Raissi, P. Perdikaris, and G. Karniadakis. Physics-informed neural networks: 477 A deep learning framework for solving forward and inverse problems involving 478 nonlinear partial differential equations. Journal of Computational Physics, 479 378:686–707, Feb. 2019. 480 48. J. C. Rajapakse, F. Kruggel, J. M. Maisog, and D. Yves von Cramon. Modeling 481 haemodynamic response for analysis of functional mri time-series. Human Brain 482 Mapping, 6(4):283–300, Jan. 1998. 483 49. W. Richter and M. Richter. The shape of the fmri bold response in children and 484 adults changes systematically with age. NeuroImage, 20(2):1122–1131, Oct. 2003. 485 50. P. N. Rosa, P. Figueiredo, and C. J. Silvestre. On the distinguishability of hrf 486 models in fmri. Frontiers in Computational Neuroscience, 9, May 2015. 487 51. M. Sarabian, H. Babaee, and K. Laksari. Physics-informed neural networks for 488 brain haemodynamic predictions using medical imaging. IEEE Transactions on 489 Medical Imaging, 41(9):2285–2303, Sept. 2022. 490 52. Z. Y. Shan, M. J. Wright, P. M. Thompson, K. L. McMahon, G. G. A. M. 491 Blokland, G. I. de Zubicaray, N. G. Martin, A. A. E. Vinkhuyzen, and D. C. 492 Reutens. Modeling of the haemodynamic responses in block design fMRI studies. 493 J. Cereb. Blood Flow Metab., 34(2):316–324, Feb. 2014.494 53. P. Sharma, W. T. Chung, B. Akoush, and M. Ihme. A Review of495 Physics-Informed Machine Learning in Fluid Mechanics. Energies, 16(5):21, Jan. 496 2023. Number: 5 Publisher: Multidisciplinary Digital Publishing Institute. 497 54. K. E. Stephan, N. Weiskopf, P. M. Drysdale, P. A. Robinson, and K. J. Friston. 498 Comparing haemodynamic models with dcm. NeuroImage, 38:387–401, 11 2007. 499 55. H. Tarbiyati and B. Nemati. Weight initialization algorithm for physics-informed 500 neural networks using finite differences | Engineering with Computers. 501 Engineering with Computers, 2024. 502 56. A. Urso, A. Fiannaca, M. La Rosa, V. Rav` ı, and R. Rizzo. Data mining: 503 Classification and prediction. In S. Ranganathan, M. Gribskov, K. Nakai, and 504 C. Sch¨ onbach, editors,Encyclopedia of Bioinformatics and Computational 505 Biology, pages 384–402. Academic Press, Oxford, 2019. 506 57. X. Wang, L. Zhou, M. He, Z. Peng, and J. Zhang. Xgboost with multi-feature 507 fusion for hemodynamic state prediction. Neuroscience, 591:63–75, 2025. 508 18/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint 58. C. L. Wight and J. Zhao. Solving allen-cahn and cahn-hilliard equations using 509 the adaptive physics informed neural networks. Communications in 510 Computational Physics, 29:930–954, 3 2021. 511 59. C. Zayane and T. M. Laleg-Kirati. A sensitivity analysis of fmri balloon model. 512 Computational and Mathematical Methods in Medicine, 2015, 2015. 513 60. S. Zeng, Z. Zhang, and Q. Zou. Adaptive deep neural networks methods for 514 high-dimensional partial differential equations. Journal of Computational Physics, 515 463, 8 2022. 516 61. S. Zeng, Z. Zhang, and Q. Zou. Adaptive deep neural networks methods for 517 high-dimensional partial differential equations. Journal of Computational Physics, 518 463:111232, Aug. 2022. 519 Supporting Information 520 Additional supporting information can be found online in the Supporting Information 521 section. 522 Balloon model constants 523 Variable Value λfin 0.2 λm 0.05 κ 1/1.54 γ 1/2.46 ‘dt‘(time step) 0.01 τM T T 30 τm 20 α 0.4 V0 0.03 ϑ0 40.3 E0 0.32 TE 0.04 ϵ 1.43 r0 25 Table 4.Model variables and typical values Residuals and Time standardisation 524 Before computational implementation, it is important to note that the orders of 525 magnitude of the residuals from equations 2 and 1 (e.g. fin, m, v, q ∼ 10−2) differ 526 significantly from the difference between data and 7 (e.g. ˆY ∼ 100); and from the 527 residuals against the initial conditions (e.g. 528 fin(t < t I); m(t < t I); v(t < t I), q(t < t I) ∼ 10−6). During gradient descent 529 optimisation, this disparity introduces complications, as it directly modulates the 530 magnitude of gradients for the loss function Ltot during backpropagation. 531 Thus, we add a scalar ponderation to level both Leq and Lic to the magnitude of 532 Ldata [25]. Before training, we standardise the input features by mapping the HRF’s 533 19/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint physical time domain [1,30] to the [−1.73,1.73] interval with unit variance, a 534 well-established technique for improving neural network convergence. Thus 535 ˚t : [0, 30] − →[−1.73, 1.73] (9) t − →t − ¯t 30·σ t (10) Where ¯t, σt are the mean and standard deviation of the original temporal coordinate. 536 Then the change of Variables Theorem compels us to rewrite the differential equations 537 accordingly 538 dq d˚t d˚t dt = 1 τM T T ·[m(˚t)− q(˚t) v(˚t) ·f out(v,˚t)] dv d˚t d˚t dt = 1 τM T T ·[f in(˚t)−f out(v,˚t)] fout(v) = v 1 α (˚t) +τ dv d˚t d˚t dt d2fin d˚t2  d˚t dt 2 +κ· d f d˚t d˚t dt =λ f ·I( ˚t)−γ · (fin(˚t) − 1) d2m d˚t2  d˚t dt 2 +κ· dm d˚t d˚t dt =λ m ·I( ˚t)−γ · (m(˚t) − 1) (11) 20/20 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716499doi: bioRxiv preprint