KurisutinaDavide Momi, Zheng Wang and John D. Griffiths. eLife 12, e83232. Version of Record published 21 April 2023. DOI; publisher PDF; PMC article and published review exchange. The 2022 preprint date is not the publication year used here.
Read on 29 September 2026. Full reading complete: all 33 pages of the final PDF, including Methods, equations, references, Appendix 1 and Appendix 2; all six main figures and ten appendix figures visually; the entire fitted-parameter CSV and five-page MDAR checklist, including its table layouts visually. Also read the entire published editor evaluation, decision letters and author response, including both decision-letter images, all 19 author-response images and the response table. The latter contains scientifically consequential clarifications absent from, or inconsistent with, the final PDF. The complete package was acquired before substantive reading. No raw-data reanalysis, independent simulation or repository code audit was performed; embedded code screenshots in the response were inspected as published evidence, not executed. References were read as bibliography, not followed as additional evidence. Archive and exact provenance.
Question and strongest supported result
Can a cortical neural-mass network reproduce individual TMS-evoked EEG potentials, and do later components in the fitted network depend on communication with other cortical regions? The authors fit a dynamical model to each person's averaged response and then remove selected connections in that model.
The result is a useful computational demonstration: an anatomically constrained, recurrent cortical model can reproduce several features of the measured response, and its later activity decreases after early disconnection of highly activated nodes. It is not a prospective prediction of an unseen perturbation, an independently recovered individual connectome, or a biological lesion experiment. No conscious versus unconscious comparison is performed.
People, anatomy and observation model
The principal TMS-EEG sample contains 20 healthy participants, mean age 24.50 years, SD 4.86, with 14 women. Existing data stimulated left primary motor cortex. A separate control dataset supplies six subjects with high- and low-amplitude motor-evoked-potential conditions. Those responses are separately fitted, rather than predicted using parameters frozen from the first cohort. The original acquisition papers were not independently read here; acquisition claims are limited to this paper's description.
The structural prior comes from 400 different HCP participants, ages 21–35, including 170 males. Their diffusion tractography uses 10 million streamlines per participant, anatomical constraints and SIFT2 weighting. The authors average the resulting weighted connectivity and tract-length matrices into a population prior, using 200 Schaefer cortical parcels and seven Yeo networks. This is MRI-scale interregional connectivity, not synaptic wiring. There is no same-person tractographic measurement for the TMS subjects in this analysis.
The connectivity matrix is transformed using its Laplacian and scalar normalization; the response specifies the Frobenius norm. Individual fits adjust connection gains around this group prior. Crucially, the author response corrects the interpretation of the stated variance 1/50: it applies to a gain in
SC_updated = SC_prior × exp(edge_gain)
(elementwise), not directly to the small raw connectivity weights. This preserves signs and zeros of the corresponding prior entries while allowing positive weights to increase or decrease. The authors' calculation of roughly +32% at two prior standard deviations is an illustrative prior scale, not a hard maximum: a Gaussian prior is not bounded at two standard deviations. Visual preservation of broad matrix patterns does not establish recovery of an individual's anatomy.
Head geometry is also largely templated: source reconstruction uses fsaverage, a three-shell boundary-element forward model and dSPM; TMS field modeling uses MNI152. The field is averaged over parcels and thresholded at 83% of its maximum. The response identifies five stimulated parcels: three left somatomotor parcels, one left dorsal-attention parcel and one left salience/ventral-attention parietal-opercular parcel. Appendix 1 specifies five fixed tissue conductivities and a 10-mm coil-to-cortex distance. Its reference to “our MRI images” should not override the explicit standard-space specification.
Thus apparent individualization comes predominantly from fitting EEG. Adjustable edge gains, two synaptic-rate parameters, four local coupling gains, global coupling, mean input and EEG magnitude scaling can accommodate response differences that might arise from physiology, anatomy, stimulation geometry or measurement. They do not uniquely assign those differences to the correct source.
Dynamical assumptions and fitting
Each of 200 cortical nodes contains a Jansen–Rit neural mass with pyramidal, excitatory-interneuron and inhibitory-interneuron populations: six first-order state equations, with local positive and negative feedback. The network has delayed interregional coupling, described as approximately 5–50 ms and derived from tract length and a conduction-velocity parameter. The final model contains no explicit thalamic or other subcortical nodes. The response describes an exploratory larger model that was not retained; this does not establish that subcortical contributions are biologically dispensable.
TMS is injected into the excitatory population. The response specifies a square input at time zero lasting 10 ms, interpreted as a coarse population charging interval rather than a literal representation of a roughly 200-microsecond pulse. EEG is derived from the difference of excitatory and inhibitory population potentials and a lead-field projection. Simulated EEG is then source reconstructed for source-space comparisons.
For each person, the algorithm repeatedly traverses one 64-channel, trial-averaged 400-ms response, from −100 to +300 ms, in sequential nonoverlapping 20-ms batches. A 20-ms burn-in precedes the baseline. PyTorch automatic differentiation and ADAM optimize parameters until convergence; the final simulation uses parameter values averaged over the last 100 batches. The Methods define mean squared error plus Gaussian parameter regularization, with distribution hyperparameters fitted. Figure 6 instead calls the loss cosine similarity, while one response passage says Pearson correlation was the fitting objective. The implemented objective cannot be resolved conclusively from these inconsistent descriptions alone.
The CSV has 20 rows and nine parameter groups: a, b, c1–c4, g, k, and mu, each with prior mean/variance and fitted mean/variance. The response defines k as EEG magnitude scaling and mu as mean input firing rate. It defines a and b as reciprocal time constants; Appendix 2 Figure 3 explicitly labels them a=1/τe and b=1/τi. Consequently, the paper's loose description of increasing b as increasing a “time constant” must not be translated directly into longer synaptic decay. Priors derive from conventional Jansen–Rit settings; the variance choices came from exploratory model development, not independent physiological measurement of these people.
Fit quality versus prediction
Figure 2 and Appendix 2 Figure 1 compare fitted and observed sensor waveforms, including all 20 individuals. Figure 3 shows grand-mean and individual source and network comparisons. The paper reports significant channel fits, significant correlations at 75.63% of source vertices, and network-level R² around .38–.46. These comparisons support reproduction of features of the fitted responses; absolute amplitudes and earliest components are visibly less well matched in some examples.
The 1,000 time-wise permutations are significance/null comparisons. They are not held-out trials or windows. Shuffling time also changes temporal structure, so rejection of that null is not evidence that the fitted model predicts an unseen physiological condition. Projecting the same fitted EEG through a source inverse supplies another representation of its agreement, not an independent source measurement. Cosine similarity, Pearson correlation and R² are used in different panels and passages and should not be pooled as one accuracy metric.
Appendix 2 Figures 6–9 extend fitting and the virtual-lesion pattern to six new subjects, separately using high- and low-motor-response data. This is a useful replication of the modeling procedure on another cohort. It does not show that parameters learned in one person, trial set, stimulation site or motor-response condition predict the other without refitting. The low-motor-response and improved sensory-masking controls address important peripheral/auditory confounds, but do not prove that every fitted late component is free of sensory contributions. Only M1 stimulation is analyzed.
A response-only stimulus-off example uses fitted parameters to generate ongoing activity. Its spectrum does not recover the empirical approximately 10-Hz peak, which the authors acknowledge. This is evidence of the present model's limited transfer beyond its fitting target, not validation of a complete resting-state model.
What the virtual lesions establish
Figures 1 and 4 explain and display the disconnection experiment. At nominal times 20, 50 and 100 ms, the authors select maximally activated nodes and set their incoming and outgoing connection weights to zero, retaining other fitted parameters. The text describes selection using the top 1% of nodes exceeding two SD above regional activity; its temporal-window specification is not entirely clear. The response reports 92% overlap with the field-defined stimulated set and says lesions act where delayed incoming signals are collected, so previously emitted signals can be blocked on arrival.
Early disconnections reduce later activity around the approximately 100-ms response and its spatial spread. Earlier/local components and some local later activity remain; the result should not be rewritten as an exact universal time dividing purely local from purely recurrent activity. Main-model network and lesion-time effects are large, including a reported interaction F(18,342)=23.79, p<.0001, partial η²=.55. These are statistics over simulations fitted to subjects, not over independently delivered biological lesions. The response adds one-subject lesions at 10-ms increments through 200 ms, supporting the qualitative timing pattern within that model.
This establishes a causal dependence of the fitted simulation on selected interregional connections under its assumed equations and parameters. It does not selectively remove feedback while preserving feedforward input, identify a unique physiological explanation, or demonstrate that a competing local/subcortical model could not reproduce the same data after fitting. There is no post-lesion refit or independently measured response to the virtual intervention. Cutting both input and output can also suppress propagation by construction. Anatomical constraints and successful fitting make the mechanism plausible; they do not convert that model intervention into a biological causal test.
Physiological parameters and identifiability
Figure 5 and Appendix 2 Figure 5 characterize two dominant spatial modes. Together they explain 74.14% of simulated and 66.96% of empirical variance. Reported peak timings are approximately 70/117 ms in simulation and 72/115 ms empirically. Across subjects, the inhibitory parameter relates negatively to the first-mode amplitude (R²=.27, p=.02) and positively to the second (R²=.28, p=.02). Neither survives the stated Bonferroni threshold .007. These are exploratory associations, despite stronger language in the abstract and discussion. Changing that parameter in a fitted simulation illustrates sensitivity, not independent identification of a person's inhibitory physiology.
Appendix 2 Figure 3 gives fitted parameter distributions; Figure 4 gives all 20 fitted matrices. The full response adds 100 repeated fits for one example person. A reviewer notes that some within-person ranges approach the between-person ranges, e.g. a approximately 94–102 across reruns versus roughly 92–102 across people. The authors argue that the distributions are concentrated and the variations small relative to parameter magnitudes, and state that final parameters average the last 100 batches. These are the authors' reported results, not independently reproduced here. Averaging batches within one fit is not equivalent to averaging 100 independent fits. Neither narrow-looking distributions nor holding parameters fixed during an ablation establishes parameter identifiability; uncertainty in mechanism conclusions across alternative equally good fits remains insufficiently characterized.
Appendix 2 Figure 10 reports fitted-connectome modularity versus response engagement, R²=.52, p=.02. The response and its code images clarify that the response quantity is derived from model-generated activity. This is a relationship among fitted model properties, not independent confirmation against the TMS participants' measured structural connectivity. No parameter-recovery experiment or anatomy-ground-truth test appears in the complete package.
PCI and consciousness
Figure 2 reports a strong association between empirical and simulated values labeled PCI, R²=.80 (text p<.0001; figure p<.001). The Methods call the metric PCI, not PCIst, describe Lempel–Ziv complexity and refer the implementation to Casali et al. (2013), without fully specifying the binarization, normalization, thresholds or calibration used here. Displayed PCI values are on a scale of tens. They should not be assumed directly interchangeable with a clinically calibrated normalized PCI threshold without resolving implementation and scaling.
The published exchange explicitly acknowledges that different waveforms can have the same PCI, and treats it as one complementary fit statistic. Agreement on this scalar is not full spatiotemporal equivalence. The study contains no unconscious-state comparison, experience-report labeling, clinical diagnostic threshold validation, or test that the simulation itself has experience. R²=.80 is not “80% of consciousness explained.”
Reproduction caveats visible in the complete package
Several descriptions remain inconsistent despite the response's clarifications: MSE versus cosine/Pearson fitting loss; raw-weight versus log-gain prior variance in the main narrative; reciprocal rates called time constants; and lesion duration described as a window in Methods but the simulation duration in one caption. Equation 4 in the final PDF visibly prints 1 − exp(...) in the denominator of a function called sigmoid. That expression is singular at the stated half-maximum voltage; the reviewer identified the sign error and the authors said it was corrected, but it persists in the archived PDF. These observations concern reporting; they are not proof that the implementation used the printed error. Resolving them would require a version-specific code audit.
The MDAR checklist mainly points to Methods/data/code sections and adds no independent validation. The paper links PyTepFit and the archived revision swh:1:rev:4222cd27fc3a451a9b43eb788d5f5e50312aed41. The exact Software Heritage directory/snapshot URL and response code-line links are preserved in provenance. These code sources and the raw datasets were not downloaded or executed for this reading; the moving main branch should not be assumed identical to that archived revision.
Implications for an individual brain model
Our inference: this paper supports combining a structural prior with explicit dynamics and perturbational observations to constrain a model. It also shows why fitting one evoked signature leaves substantial ambiguity: estimated wiring, local physiology, inputs and observation scaling can compensate for one another. The study does not test anatomy alone, establish unique recovery of hidden parameters, or show that the fitted system generalizes across states and interventions.
A stronger individual-model test would freeze parameters after calibration, predict new stimulation sites, intensities and sessions, compare same-person anatomy with group priors, and test whether equally good initial fits diverge on those held-out outcomes. It would propagate fitting uncertainty through the lesion conclusions and separately validate any consciousness-related metric. These are proposed next tests, not results of this paper. Nothing here establishes preservation of autobiographical memory, personal identity, subjective experience or a conscious emulator.
Code audit (added 2026-09-29, main session)
A read-only audit of PyTepFit at the paper's archived revision 4222cd27 settles several caveats above: report. The git tree hash matches the cited Software Heritage directory.
- PCI. The reported "PCI" is PCIst (Comolatti et al. 2019), not Casali's Lempel–Ziv PCI. For the simulated values, the simulated response is pasted over the empirical grand-average baseline. The released cell labels Pearson r as "R2".
- Fitting. The loss is RMSE. The code's sigmoid is the correct logistic, so Eq. 4's
1 − expis a typesetting error. - Priors. The connectivity prior is absent from the loss. The parameter prior's fitted mean tracks the parameter, so the prior does not constrain it. The leadfield is fitted.
- Timing. Batches are 50 ms, not 20 ms. The stimulus sits at +10…+20 ms.
- Supplementary file 1. It reports last-batch values and a fitted precision hyperparameter, not posterior summaries.
- Evaluation. All fit statistics are in-sample.
- Reproducibility. The release cannot reproduce the paper as shipped: inputs are missing and paths are hard-coded.