- 1Advanced Institute for Materials Research (WPI-AIMR), Tohoku University, Sendai, Japan
- 2Departament de Física de la Matèria Condensada, Universitat de Barcelona, Barcelona, Spain
- 3Universitat de Barcelona Institute of Complex Systems (UBICS), Barcelona, Spain
- 4Research Institute of Electrical Communication (RIEC), Tohoku University, Sendai, Japan
- 5Faculty of Science and Technology, Oita University, Oita, Japan
- 6The School of Systems Information Science, Future University Hakodate, Hakodate, Japan
- 7International Research Center for Neurointelligence (WPI-IRCN), The University of Tokyo, Tokyo, Japan
Impaired brain function is restored following injury through dynamic processes that involve synaptic plasticity. This restoration is supported by the brain’s inherent modular organization, which promotes functional separation and redundancy. However, it remains unclear how modular structure interacts with synaptic plasticity to define damage response and recovery efficiency. In this work, we numerically modeled the response and recovery to damage of a neuronal network in vitro bearing a modular structure. The simulations aimed at capturing experimental observations in cultured neurons with modular traits which were physically disconnected through a focal lesion. The damage reduced the frequency of spontaneous collective activity events in the cultures, which recovered to pre-damage levels within 24 h. We rationalized this recovery in the numerical simulations by considering a plasticity mechanism based on spike-timing-dependent plasticity, a form of synaptic plasticity that modifies synaptic strength based on the relative timing of pre- and postsynaptic spikes. We observed that the in silico numerical model effectively captured the decline and subsequent recovery of spontaneous activity following the injury. The model supports that the combination of modularity and plasticity confers robustness to the damaged neuronal network by preventing the total loss of spontaneous network-wide activity and facilitating recovery. Additionally, by using our model within the reservoir computing framework, we show that information representation in the neuronal network improves with the recovery of network-wide activity.
1 Introduction
Brain injury impairs crucial brain functions, including cognition, memory, and higher-level executive functions (Walker and Tesco, 2013). Although injuries may temporarily affect large areas of the brain, these functions are often restored spontaneously or through rehabilitation within a few months (Chen et al., 2010; Christensen et al., 2008; Edlow et al., 2021; Murata et al., 2015; Nudo, 2013). Restoration of neuronal activity is central to the recovery of brain function, a process that has been extensively investigated in vivo (Fukui et al., 2020) and in vitro (Ayasreh et al., 2022; Teller et al., 2020). This restoration is believed to be mediated by network-wide reorganization of functional connections through synaptic plasticity. Clinical and computational studies have indicated that spike-timing-dependent plasticity (STDP), a form of Hebbian synaptic plasticity, is a potential mechanism for restoring damaged neuronal networks (Bunday and Perez, 2012; Gabrieli et al., 2020; O’Neill et al., 2023). Taken together, these findings suggest that synaptic plasticity restores neuronal activity, which in turn leads to the recovery of higher brain function.
In addition to plasticity, inherent network topological traits are highly involved in injury response and recovery (Arnemann et al., 2015; Boroda et al., 2021). Recent brain connectome analyses have revealed that the brain network evolutionarily conserved a modular structure in which densely connected neuronal populations (modules) are relatively sparsely connected to other modules (Lee et al., 2016; Lynn and Bassett, 2019; Meunier et al., 2010). The coexistence of highly clustered modules and shortcuts between modules facilitates redundancy, ensuring resilience to failure by diversifying information flow, thereby granting robustness to network functions. For example, animal studies have shown that modularity of the frontal cortex increases task feasibility under partial perturbations (Chen et al., 2021). Furthermore, modularity is advantageous for functional recovery, as patients with higher modularity of the cortex exhibit greater improvements in executive function during cognitive training after brain injury (Arnemann et al., 2015). However, it remains unclear how the underlying modular structure of the neuronal network combines synaptic plasticity, prominently STDP, to lay out efficient damage response and recovery.
In this study, we conducted experiments on cultured neurons and a spiking neural network (SNN) model, both of which were designed with a modular organization, to examine their response to damage. We created modular neuronal cultures using topographically modulated substrates shaped as parallel tracks, which effectively guided and constrained connectivity along tracks (Montalà-Flaquer et al., 2022) and later inflicted focal damage using a scalpel. Immediately following injury, the rate of spontaneous neuronal activity was substantially reduced but recovered to its original level after 24 h. Next, an SNN model was devised, and its parameters adjusted to align with experimental observations, resulting in a similar reduction and recovery of network activity after damage. The constructed SNN was then used to investigate the effects of damage at different module locations, and to examine the response to damage in networks without a modular structure. The in silico replication of the damaged networks suggested that the underlying modular organization helps preserve overall network organization upon damage, resulting in faster recovery compared to non-modular networks. Finally, the recovery of neuronal activity was conceptually linked to the recovery of brain function using a spoken digits classification task within a reservoir computing framework. Our results underscore the interplay between network architecture and plasticity in enabling the restoration of activity and basic functionality following damage, providing new insights into recovery mechanisms in living neuronal assemblies.
2 Materials and methods
2.1 Fabrication of engineered topographical substrates
Polydimethylsiloxane (PDMS) topographical substrates (Montalà-Flaquer et al., 2022) were used to fabricate cultured neuronal networks with modular characteristics. The substrates were prepared using a specially designed mold made of fiberglass and copper (2CI Circuitos Impresos, Spain) to form a two-layer structure. The bottom layer consisted of a uniform fiberglass sheet 2 mm thick, while the top layer contained copper motifs 70 μm high, shaped as parallel stripes 300 μm wide, separated by 300 μm fiberglass gaps, and extending the entire length of the mold. In the remainder of this paper, we refer to this pattern as “tracks.”
PDMS (Sylgard 184, Ellsworth Adhesives) was poured onto the fiberglass-copper mold as a mixture of 90% base and 10% curing agent, and cured at 90°C for 2 h. The PDMS cast was then gently peeled off to form a two-level substrate, which was the negative of the original mold, with PDMS valleys corresponding to copper imprints and crevices to the fiberglass, respectively. The PDMS sheet was perforated using a stainless-steel punch (Bahco 400.003.020), resulting in disks 6 mm in diameter and 1 mm in thickness. The topographical “tracks” pattern on the surface contained 300 μm wide modulations extending across the entire disk. Conceptually, each track gave rise to a module in which neurons were strongly connected along the tracks and weakly connected across them. Before culturing, the PDMS disks were cleaned with ethanol, dried, and mounted on clean coverslips (#1 Marienfeld Superior; 13 mm diameter). One or two PDMS disks were placed on each coverslip, and the PDMS-glass assembly was sterilized in an autoclave (Selecta 4002515). Sterilization increased the bonding between the PDMS and the glass surface such that the PDMS remained attached to the glass throughout the lifespan of the culture.
2.2 Cell culture and GCaMP6s viral transduction
Primary neurons derived from rat embryonic cortices on days 18–19 were used in all experiments. Animal experiments and tissue manipulations were conducted following the approval order B-RP-094/15-7125 from the Ethics Committee for Animal Experimentation of the University of Barcelona (CEEA-UB). Rats were provided by the Animal Farm of the University of Barcelona. Brain dissection was performed in ice-cold L-15 medium (Gibco), and cortical tissue was mechanically dissociated by pipetting. Neurons were then transferred to a “plating medium” [90% Eagle’s minimum essential medium (MEM, Invitrogen), 5% horse serum (HS, Invitrogen), 5% bovine calf serum (Invitrogen), 1 μL/mL B27 (Sigma)]. Before seeding the neurons, the PDMS topographical surfaces were coated overnight with a solution of 20 μg/mL poly-L-lysine (PLL; Sigma-Aldrich) in borate buffer. Neurons were seeded onto these surfaces to create cultured neuronal networks at a density of approximately 400 neurons/mm2. On day in vitro (DIV) 1, neurons were transduced with the genetically encoded fluorescence calcium indicator GCaMP6s (AAV9.Syn.GCaMP6s.WPRE.SV40, Addgene) under the Synapsin 1 promoter, so that only mature neurons (and not other cells such as glia) expressed the fluorescence indicator. On DIV 5, the “plating medium” was replaced with a “changing medium” (90% MEM, 10% HS, 0.5% 5-Fluoro-2-deoxyuridine) to limit glial growth. On DIV 8, the medium was switched to “final medium” (90% MEM and 10% HS), which was periodically refreshed every 3 days. Cultured cells were incubated at 37°C, 5% CO2, and 95% humidity (Memert INCO2-246).
2.3 Monitoring of neuronal activity
Wide-field calcium imaging of the PDMS topographical cultures was performed using an inverted fluorescence microscope (Zeiss Axiovert 25C, Zeiss GmbH) equipped with a high-speed camera (Hamamatsu Orca Flash 4.0v3, Hamamatsu Photonics) and a fluorescent light source (mercury vapor arc lamp, Osram GmbH) at 12–13 DIV. The fluorescence image series was recorded at 50 frames/s, 8-bit grayscale levels, with an image size of 1,024 × 1,024 pixels. The recordings covered a field of view of 7.1 × 7.1 mm2, which was achieved by combining a 2.5× objective with an optical zoom, allowing observation of the entire 6 mm diameter culture. Recordings of spontaneous activity were 15 min long and were controlled through the Hokawo 2.10 software (Hamamatsu Photonics). During recordings, the cultured neuronal networks were placed in a glass microincubator (Ibidi GmbH), which maintained the same environmental conditions as the standard incubator. The temperature was set to 25°C to facilitate spontaneous activity.
2.4 Induction of damage in the neuronal cultures
The focal injury was induced through a scalpel on DIV 12−13 by inflicting a 4-mm long straight incision on the neuronal population in the 6 mm diameter PDMS. The protocol for monitoring network behavior through a complete experimental “damage sequence” was as follows: spontaneous activity was recorded for 15 min before injury. Immediately after injury, neuronal activity was recorded for 30 min; and thereafter for 15 min at quasi-logarithmic time points post-injury, specifically at 2, 6, and 24 h post-injury (Ayasreh et al., 2022).
2.5 Activity quantification
On the images of acquired recordings, a series of 150 × 150 μm2 square regions of calcium intensity data were defined as regions of interest (ROIs), resulting in a total of 1,400 regions. Spike trains of neuronal activity were extracted from each region by first generating a time series of mean fluorescence intensity within each ROI and then by detecting inferred spiking events using a Schmitt-trigger filter (lower threshold = 0.4 and upper threshold = 0.8) (Soriano, 2023). ROIs with spike rates below 0.01 Hz were considered to contain no active cells and were excluded from the analysis.
Neuronal activity comprised either individual activation of ROIs or collective events in which several ROIs displayed coordinated activity within a short time window. These coordinated events were termed network bursts and reflect the synchronized and rhythmic activity of neuronal populations in vitro (Wagenaar et al., 2006; Yamamoto et al., 2016; Pasquale et al., 2017). In our study, network bursts were considered significant when the number of participating ROIS within a 200 ms window exceeded 20% of the total number of ROIs, as below this fraction a burst could not be easily distinguished from random background activity. The number of bursts during a recording and their amplitudes (network fraction) were used to quantify neuronal culture activity before and after injury, as well as during recovery.
2.6 Effective connectivity
Causal interactions between pairs of active ROIs were inferred using generalized transfer entropy (TE) (Stetter et al., 2012; Tibau et al., 2020; Montalà-Flaquer et al., 2022). Given the spike trains of a pair of ROIs X and Y, with X = xm and Y = ym, the amount of information transferred from ROI Y to X is given by:
where m is a discrete-time index and k ( = 2) is the Markov order. For inference, the 15-min-long spike train of each ROI was built using a time bin of 20 ms, with data containing either “1” for the presence of a spike or “0” for its absence. Transfer entropy is a nonlinear and nonsymmetric measure in X and Y (TEY→X≠TEX→Y), allowing us to estimate causal interactions, i.e., effective connectivity, in the network. Significance was established by first normalizing the distribution of the TE values using a z-score transformation:
where ⟨TE⟩ is the mean value of all TE scores and σ their standard deviation, and then by setting a threshold zth = 2 so that zY→X = 1 ∀ zY→X≥zth, and zY→X = 0 otherwise. The final adjacency matrix Z of the effective connections was directed and binary.
2.7 In silico neuron model
The Izhikevich model was used in this study (Izhikevich, 2003) since it accurately reproduces the spiking behavior of cultured neurons with sufficient biological plausibility (Orlandi et al., 2013). The model is described by the following equations:
where v is the membrane potential, u is the recovery variable, IE is the excitatory synaptic input, II is the inhibitory synaptic input, and Iin the external current. During spontaneous activity, an external current was absent, indicating that Iin = 0. The membrane potential was reset when it exceeded vpeak = 30 mV, at which point an action potential (a spike) was recorded for that neuron. The parameters a, b, c, and d determine the dynamic characteristics of the modeled neurons. In the present study, the simulated neuronal networks contained 80% excitatory and 20% inhibitory neurons. For excitatory neurons, the parameters were set as [a, b, c, d] = [0.02, 0.2, −65.0, 8.0], whereas for inhibitory neurons they were set as [a, b, c, d] = [0.1, 0.2, −65.0, 2.0].
2.8 In silico network generation and infliction of damage
Neurons were randomly positioned in a homogeneous manner within a circular area of 3 mm in diameter, with a density of 400 neurons/mm2, resulting in approximately 2,800 neurons. Following the method of Orlandi et al. (2013), dendrites were modeled as circles centered on each neuronal soma with a radius of 150 μm, while axons were described through a growth process as follows. For each axon, its maximum length ℓ was first determined by sampling from a Rayleigh distribution with an average length of ⟨ ℓ ⟩ = 1.1 mm. Then, segments of 0.1 mm in length were placed in a concatenated manner along a pseudo-straight path, such that, at each growing step, a new axon segment slightly deviated from the direction of the previous one by an angle θi, according to the following probability:
where i represents each growing step, and σθ was set to 5.73° (0.1 rad). The effect of PDMS on the development of the neuronal culture was modeled following (Houben et al., 2025). To account for the PDMS modulations, present in the living neuronal networks, virtual parallel bands 0.2 mm wide were considered for the “crevices” of the PDMS, which were separated by 0.3 mm wide “valleys.” These bands acted as obstacles that interfered with the axonal growth, such that axons crossed from top to bottom with a 50% probability (Pdown in Figure 1C) and from bottom to top with a 5% probability (Pup). Whenever an axon failed to cross, it was deflected and continued to grow parallel to the obstacle wall. Simulations also included “control” scenarios with no virtual bands to investigate the impact of spatial constraints on network dynamics.

Figure 1. Schematic pipeline of the experimental and numerical approaches. (A) Primary cortical mouse neurons were cultured on an engineered substrate shaped with PDMS topographical modulations. The topographical substrate was a 6-mm-diameter, 1-mm-thick disk. Its surface featured 70 μm-high parallel stripes, each 300 μm wide and spaced 300 μm apart. A bright-field image of a prepared neuronal culture is shown on the right. Network structure under the same culture conditions was previously validated by immunostaining (see Montalà-Flaquer et al., 2022). (B) Damage was applied to the culture with a scalpel, creating a wound approximately 4 mm long (blue box) at day in vitro 12–13. Spontaneous activity was then recorded using calcium imaging just before the damage and at preset time points: immediately after the damage and at 15 min, 2 h, 6 h, and 24 h post-injury, to monitor changes in activity levels and other properties. (C) In the in silico modeling, an axon growth model was used to replicate the structure of the cultured neuronal networks. Pup and Pdown represent the probability that an axon crosses the patterned bands upwards or downwards, respectively, and are set to 5 and 50%. Then, a spiking neural network model was initialized, and synaptic weights were adjusted using STDP over 72 h. Injury was applied in the same manner as in the experiments, and network activity was evaluated by simulating STDP again at the same time points as in the cultured experiments. Scale bars in (A,B) are 1 mm. Cartoons in (A,B) were created with BioRender.com.
In both the “control” and “tracks” scenarios, a connection between neurons i→j was established whenever the axon of neuron i crossed the area covered by the dendritic tree of neuron j with a 20% probability. The resulting binary structural connectivity matrix A = {aij} (1 for the presence of a connection and 0 otherwise) was then used to generate a new matrix W = {wij} that preserved the connectivity relationships but incorporated weights. The specific values of wij are described in section 2.10.
To model damage in the in silico networks, a straight line was defined as a reference in the center of the simulated culture, placed horizontally with a length of 1.5 mm (half the diameter of the network). Subsequently, all connections crossing this line were set to 0 in the connectivity matrix A. Additionally, neurons whose axons were severed due to damage were considered dead (to mimic axonal transection) and excluded from further membrane potential updates during the simulation of the dynamics.
2.9 Synaptic model
The numerical simulations utilized the following synaptic model, which is characterized by an exponential decay of induced postsynaptic membrane currents:
In this model, j is the postsynaptic neuron index, the time of the k-th somatic spike of presynaptic neuron i, and δ(⋅) the Dirac delta function. The set Exc refers to the excitatory neurons, and Inh refers to the inhibitory neurons. The average excitatory weight, can be modified to alter the capacity of the network to display synchronized activity events (network bursts). The parameters τE and τI account for the time constants of current loss and were set to τE = 5 ms for excitatory synapses and τI = 20 ms for inhibitory ones. Additionally, axonal conduction delays dij on excitatory connections were uniformly distributed in the range of [0, 5] ms, while those on inhibitory connections were fixed to 1 ms. The total excitatory and inhibitory input onto a neuron, IE and II, respectively, was the sum of synaptic currents from incoming input connections. Additionally, the target neuron j receives Poisson noise in the form of spontaneous synaptic inputs at an average rate of 1 Hz. The quantity tjp denotes the time of the p-th Poisson spike input to neuron j, and ξ is the noise amplitude. The higher the frequency or strength of the Poissonian noise, the higher the capacity of the network to spontaneously activate, either at a single neuron level or as collective network burst.
2.10 Spike-timing-dependent plasticity
The STDP model (van Rossum et al., 2000; Rubin et al., 2001; Gütig et al., 2003) is based on a phenomenon in which the strength of synaptic weight changes depending on the time difference between the firing of pre-synaptic and post-synaptic neurons. This is described mathematically by the following set of equations:
where wij is the synaptic weight from neuron i to j and Δ t = tj − ti − dij the time difference between the firing of pre-synaptic neuron i and post-synaptic neuron j. τ is a constant set to 20 ms that accounts for the characteristic time for the strengthening or weakening of a connection. A+(wij) and A−(wij) are functions of wij (van Rossum et al., 2000; Rubin et al., 2001; Gütig et al., 2003) given by:
where η+ and η– are learning constants, set to 0.1 and −0.12, respectively. The wmax is the maximum value of the synaptic weight and is set to 6.8 unless otherwise specified. With A+(wij) and A−(wij) set in this way, the average value of synaptic weights is approximately . Since only the connection strengths between excitatory neurons were updated, all other connections were set to: wEI = for excitatory-to-inhibitory connectivity, wIE = for inhibitory-to-excitatory, and wII = and for inhibitory-to-inhibitory. Nearest-neighbor spikes only contributed to modifying weights on plastic connections (Izhikevich and Desai, 2003; Morrison et al., 2008).
2.11 Computation of global efficiency for structural and effective networks
Network-wide connectivity of the in silico, weighted synaptic connectivity matrices W = {wij}, before and after injury, was quantified by means of the global efficiency Eglob (Latora and Marchiori, 2001; Onnela et al., 2005; Fagiolo, 2007; Rubinov and Sporns, 2010) as:
where N is the number of nodes in W, and lij is the shortest path length between nodes i and j, obtained from all possible path lengths Lij through Dijkstra’s algorithm. Lij was computed from the weighted connectivity matrix W as:
We note that Lij transforms weights into connection lengths with the property that stronger weights correspond to shorter path lengths, increasing global efficiency.
Similarly, for the computation of global efficiency based on the TE-derived effective connectivity matrices, the path length Lij was calculated as:
where TEi→j represents the transfer entropy values prior to z-score normalization, and TEmax is the maximum value of these unnormalized transfer entropy values observed across all time points of a given damage sequence.
2.12 Normalized mutual information
NMI is a measure of the similarity between two partitions (Alotaibi and Rhouma, 2022) and was used to quantify the differences in the community structure of the synaptic matrix W or effective connectivity matrix Z before and after damage:
where N is the total number of neurons, and C is the best partition of the matrix (either W or Z), computed through the Louvain algorithm (Blondel et al., 2008). |c| and |c’| denote the number of neurons in the communities c and c’ that belong, respectively, to partitions C and C’. |c ∩ c′| indicates the number of neurons in the intersection of c and c′. The greater the similarity in community partitions between two conditions, the higher the NMI. In the present work, NMI was used to compare the similarity in community structure between pre- and post-injury networks, with NMI = 1 indicating identical structure.
2.13 Reservoir computing
The reservoir computing numerical implementation comprised an input layer, a reservoir layer, and an output layer. Input signals were spoken digits from the TI-46 dataset (Liberman et al., 1993). The spoken digits were first converted to a 78-channel cochleagram using Lyon’s passive ear model (Slaney and Lyon, 1993). The cochleagram was then normalized to a maximum value of 1. The reservoir layer consisted of a spiking neural network, with 5% of the total neurons receiving signals from the input layer. The signals delivered to each of these neurons were generated by summing two randomly selected channels of the 78-channel cochleagram. This input was then multiplied by the averaged synaptic weight ŵ and added to Iin in Eq. (3). During the task, the updating of synaptic weights was stopped. The reservoir state x(t) was constructed from a subset of 5% of the neurons in the network, but different from the subset that received the inputs, according to the following equation:
where tik is the time of the k-th spike of neuron i. The time constant τx was set to τx = 1 s, and xstep was set to 0.1. The output y(t) was calculated as:
where Wout is the output weight matrix, which was obtained using ridge regression during the training phase as:
where X represents the reservoir state and was constructed by horizontally concatenating x(t) from the start of the first trial to the end of the last trial. Each trial consisted of an input of a spoken digit, followed by a 10-s output period. The matrix represents the target signal and was created by concatenating ŷ(t), a target vector, where the element corresponding to the correct label is set to 1 for 2.5 s following each stimulus onset and zero otherwise. The λ ( = 1) is the regularization coefficient, and I is an identity matrix. The training dataset included 10 distinct samples of each spoken digit labeled as “zero,” “one,” and “two” from the speaker identified as “f1.” The test dataset was identical to the training dataset and was used to examine the loss and recovery of function due to damage. The estimated output was obtained from y(t) as , where i is an element of the output layer, tk the onset of the k-th input, and Δt the target duration (set to 2.5 s). Accuracy was evaluated as the fraction of correct estimates.
3 Results
We investigated damage in neuronal networks and numerically explored the interplay between modularity and plasticity in promoting recovery. The results below are organized as follows. Section 3.1 provides an overview of the experimental design and describes the damage and recovery observations in both the cultured neuronal networks (in vitro) and the spiking neural network (SNN, in silico). Section 3.2 focuses on the numerical exploration of synaptic weight changes upon network damage and elucidates on the capacity of STDP to promote network reorganization and recovery in silico. Section 3.3 further utilizes the in silico model to investigate the dependence of the damage and recovery processes on the location and size of the injury, as well as the effect of modular architecture. Finally, section 3.4 integrates the in silico model into a reservoir computing framework to demonstrate that information processing capability is recovered following the recovery of the network.
3.1 Self-organized recovery after injury in cultured and spiking neuronal networks
We prepared cultured neuronal networks with a modular structure by growing rat primary neurons on topographical PDMS substrates as described previously (Montalà-Flaquer et al., 2022; Figure 1A). Crevices and valleys of the PDMS surface repeat in one direction, forming parallel tracks. There were approximately 10,000 neurons on a 6 mm diameter PDMS culture, and their behavior was monitored using calcium fluorescence imaging. The cultured neuronal networks were damaged at DIV 12 or 13 using a scalpel in a direction transverse to the PDMS tracks (Figure 1B). Subsequent recovery of neuronal activity was measured immediately after injury and at 15 min, 2 h, 6 h, and 24 h later. To mimic the experimental design, an SNN model was created to understand the synaptic and network mechanisms of the recovery process of the cultured neuronal network (Figure 1C). Network connectivity was modeled by simulating the growth process of nerve axons on a substrate that incorporated experimental-like parallel modulations, effectually shaping a modular network. The created neural network models were initialized and simulated with STDP for 72 simulated hours to settle all transients in the weight distribution. Then, the network model was damaged similarly to the culture experiments.
Figure 2A shows fluorescence images of representative cultured neuronal networks before and after damage. Neuronal activity for each recording was extracted from the average fluorescence intensity of 1,400 ROIs that covered the culture area (see Supplementary Figure S1). The cultured neuronal network exhibited spontaneous activity characterized by collective quasi-synchronous events (network bursts). Such a bursting behavior was observed both before and after damage (Figure 2B); however, network burst frequency was reduced immediately after damage. Effective connections were computed from firing patterns using transfer entropy, and functional communities were evaluated. These communities revealed groups of neurons that interacted more strongly within their groups than with the rest of the network.

Figure 2. Damage and recovery in cultured neurons and in silico. (A) Fluorescence images before (left) and after damage (right) in vitro. The blue dotted line indicates the PDMS substrate, and the red box marks the damaged area. (B) Raster plots of neuronal activity before, immediately after, and 24 h after injury in cultures. (C) Effective connectivity maps and modular partitions. Functional modules are aligned with PDMS tracks before and 24 h after damage, but become disorganized immediately after it. Module colors are not matched across time points due to post-injury reorganization. (D) NMI values between pre- and post-injury partitions in vitro (n = 4). (E) Ratio of bursting events after damage relative to before damage (n = 4). (F) Structural maps of the network before and after damage in silico. (G) Simulated raster plots before, immediately after, and 24 h after injury. (H) Effective connectivity and modular structure in silico, showing acute disruption and partial reorganization. (I) NMI values in simulations. The bars denote the mean values, and error bars indicate the 95% confidence interval (n = 15). (J) Burst activity ratio in simulations (n = 15), showing recovery similar to the experimental observations. Scale bars: 1 mm (A,C,F,H).
As shown in Figure 2C, the spatial organization of the detected communities changed substantially throughout the network before and after injury (see Supplementary Figure S2). Before damage, the functional communities were aligned along the “tracks” pattern, indicating that the substrate’s topographical modulation facilitated the formation of the communities. Such alignment was lost immediately after damage, suggesting that the local injury caused long-range alterations. The similarity of the community structure before and after injury, quantified through the normalized mutual information (NMI), decreased after damage and remained low even at 24 h post-injury (Figure 2D). Interestingly, however, the cultured neurons reorganized their functional communities differently from the pre-damage communities, recovering the distribution of connection angles after 24 h (Supplementary Figure S2B).
Since activity in neuronal cultures is characterized by network bursts, changes in the frequency of such bursting were used as the first approach to quantify the effect of damage. As shown in Figure 2E, the “burst count ratio,” that is, the ratio of the number of bursting events relative to pre-damage conditions, substantially dropped immediately after injury. Upon recovery, the ratio approached pre-damage levels or even exceeded them 24 h after injury.
To understand the synaptic mechanisms underlying this recovery, we utilized an SNN model, simulating a population consisting of 80% excitatory and 20% inhibitory neurons, where excitatory coupling was modified by STDP (Figure 2F). Raster plots of activity were then generated, and the model parameters were adjusted—specifically the synaptic weights wmax—to qualitatively reproduce the experimental observations. The variability observed in the cultured neurons was captured by modifying either wmax or the frequency of spontaneous noise (Supplementary Figure S3). In general, regardless of parameter settings, and as long as there was sufficient drive for spontaneous activity, we observed that the frequency of bursting decreased following damage and then gradually recovered. This exploration demonstrates that the recovery process is robust to variations in the parameters that control overall activity level. In the subsequent analyses, we considered numerical simulations with wmax = 6.8 and a spontaneous noise frequency of 1.0 Hz, which provides the average behavior.
Finally, we conducted the same analysis pipeline as in the experiments to investigate whether the model captured the experimental observations on activity and effective connectivity. When damage was inflicted on the SNN model, a decrease in activity was observed (Figure 2G), and, as in the culture experiment, the organization of functional communities shifted from a track-oriented to a mixed arrangement immediately after damage (Figure 2H; Supplementary Figure S4). The similarity of the community structure showed a sustained reduction post-injury (Figure 2I), whereas the distribution of connection angles recovered to a track-oriented configuration by 24 h (Supplementary Figure S4B). By qualitatively comparing with experiments, these results suggest that STDP of excitatory-to-excitatory synapses is sufficient to model the recovery of spontaneous activity and damage-induced alterations in both dynamics and functional organization. Additionally, consistent with culture experiments, a quantification of the ratio of burst events revealed a transient decrease by a factor 0.68 immediately after injury, recovering to the baseline level after approximately 24 h of simulation (Figure 2J).
We note that functional recovery in silico did not occur when STDP was absent (Supplementary Figure S5). However, such a scenario cannot be experimentally tested in cultured neurons, as plasticity is an intrinsic property of living neuronal assemblies and cannot be easily suppressed without affecting other processes important for physiological activity. Moreover, multiple plasticity mechanisms may act concurrently, e.g., synaptic scaling and homeostasis (Turrigiano, 1999, 2008; Effenberger et al., 2015).
3.2 Reorganization of synaptic weights mediates recovery of neuronal activity in silico
To investigate the potential of STDP to spatially redistribute synaptic weights and restore activity in silico, we analyzed the changes in excitatory-to-excitatory synaptic weights of the SNN model before and after injury. Figure 3A shows the change in synaptic weights of a representative network immediately after injury and 24 h later. The increase and decrease in synaptic weights occurred throughout the entire network rather than at specific locations. We observed large fluctuations in synaptic weights post-damage, whereas the undamaged network remained relatively stable (Figure 3B), suggesting that injury enhanced the effect of STDP. Despite the large variation in the damaged network, the mean synaptic weight remained practically constant with a value of about 3.1 before injury and after reorganization (Figure 3C).

Figure 3. Changes in synaptic weights in silico due to STDP. (A) Representative network illustrating the changes in synaptic weights 24 h post-damage, as Δw = w24h − winjury. Red connections represent those that increase in weight, while blue connections represent those that decrease in weight. Scalebar is 1 mm. (B) Evolution of synaptic weights changes relative to the weights immediately after injury (0 min) for the representative damaged network and control network. The plot shows the 100 representative synaptic weights at each time point, with the values calculated by subtracting the weights at 0 min. (C) Evolution of the mean value of the synaptic weights for damaged networks and control networks, showing that the global change is minimal. The blue curve corresponds to a damaged network, whereas the gray corresponds to control network. Shadings indicate the 95% confidence interval, n = 15 networks. (D) Evolution of the global efficiency of the weighted connectivity matrix W after injury and during recovery. (E) Relationship between global efficiency and bursting event rate, with data exhibiting a Pearson’s correlation coefficient of r = 0.80 and p < 0.001 [two-sided Pearson correlation test, n = 75, degrees of freedom (df) = 73].
We then investigated the effect of post-damage reorganization of synaptic weights driven by STDP on the efficiency of information transfer in a network. We observed that the global efficiency of the damaged network decayed immediately after injury and gradually recovered over time (Figure 3D, blue curve), whereas the control condition exhibited only a small variation (gray curve). Global efficiency quantifies the ease of neuronal communication across the entire network. Thus, the results indicate that communication between neurons was reduced due to damage but was restored through the reorganization of synaptic connections via STDP. Global efficiency was further examined in relation to the rate of network burst frequency, and we observed a strong correlation between the two quantities [r = 0.80, p < 0.001, two-sided Pearson correlation test, n = 75, degrees of freedom (df) = 73] (Figure 3E). This correlation indicates that synchronous activity requires reliable communication between distant neurons in the network and, therefore, high network efficiency for information exchange.
Furthermore, we quantified global efficiency of the effective connectivity in the cultured neuronal networks estimated via transfer entropy. Although the values varied across samples (Supplementary Figure S6A), we observed a consistent decrease upon damage followed by gradual recovery when each damage sequence was normalized by its pre-damage value (Supplementary Figure S6B). A similar analysis for the STDP-driven in silico model also revealed a comparable decrease and subsequent recovery (Supplementary Figure S6C). The magnitude of efficiency in silico aligns with in vitro samples with higher initial global efficiency, while those with lower efficiency may also be modeled by adjusting simulation parameters, such as connection density. These findings suggest that the network alterations observed in the experimental neuronal cultures are compatible with the numerical model that considers STDP as the main recovery mechanism. In summary, the STDP reorganization brings the post-damage network to a new state with increased global efficiency and activity level, but with different spatial distribution of connection weights.
3.3 Effect of modular organization on damage and recovery in the in silico model
We next examined the influence of the modular structure, shaped by a series of interconnected parallel tracks, to provide robustness against various damage scenarios. We designed several damage conditions, alongside control scenarios, in which either intra- or inter-modular connections were damaged, as shown in Figure 4A. Cuts perpendicular to and along the tracks affected the intra-modular and inter-modular connections, respectively. A network generated by modeling axonal growth on an unpatterned surface was used as a control. Two different damage extents were considered: one in which the length of the injury extended half the diameter of the network and another in which the injury extended the full width (see Figure 4A).

Figure 4. Importance of modularity on the magnitude and direction of damage in silico. (A) Schematic representation of damage actions on patterned (modular) and homogeneous networks. Intra-modular or inter-modular damages are achieved by, respectively, cutting the network along the direction perpendicular or parallel to the tracks’ orientation. In the sketch, “half” indicates that the cut extends half the diameter of the network, while “full” indicates that the damage effectually separates the network in two parts. “Damage on non-module” represents the control scenario of neurons simulated on a flat surface. (B) Burst frequency changes for the six damage scenarios. For each scenario, the color bar shows the ratio of burst frequency between the damaged and control conditions. Color bars use the same color scheme as in (A). The bar height indicates the mean, and error bars indicate the 95% confidence intervals. Color asterisks indicate a two-sided one-sample t-test (*p < 0.05; **p < 0.01; ***p < 0.001; N.S., no significance; n = 15, df = 14). Asterisks in black correspond to a two-sided unpaired t-test (**p < 0.01; ***p < 0.001; N.S., no significance; n = 15, df = 13). (C) Sketch illustrating the concept of normalized mutual information (NMI). Community structure is obtained via the Louvain algorithm. (D) NMI values between pre- and post-injury networks, for the six damage scenarios. Bars show mean value, and error bars indicate the 95% confidence intervals (n = 15). The horizontal line at NMI = 1 represents the reference point for identical community structures between pre- and post-injury conditions.
Figure 4B shows the rate of network burst relative to the non-damaged control at each time point after injury for a total of six conditions (four in modular networks and two in non-modular networks). Immediately after injury (0 min), the activity was lower than the control in all conditions [intra-modular half, 0.27 ± 0.14 (mean ± SD), two-sided one-sample t-test, p < 0.001 (vs. 1.00); intra-modular full, 0.06 ± 0.07, p < 0.001; inter-modular half, 0.60 ± 0.17, p < 0.001; inter-modular full, 0.41 ± 0.11, p < 0.001; non-module half, 0.27 ± 0.14, p < 0.001; non-module full, 0.09 ± 0.06, p < 0.001; n = 15, df = 14]. A comparison between half-cuts and full-cuts showed that networks with a full-cut had a significant reduction in activity [two-sided unpaired t-test, p < 0.001 (intra-modular half vs. intra-modular full); p < 0.01 (inter-modular half vs. inter-modular full); p < 0.001 (non-module half vs. non-module full)]. Additionally, damage perpendicular to the tracks direction (intra-modular) had a greater impact on the rate of neuronal activity than damage applied along the tracks (inter-modular) (two-sided unpaired t-test, p < 0.001, intra-modular half vs. inter-modular half). However, there were no significant differences between non-modular networks and modular networks in which damage was applied within a module, i.e., perpendicular to the tracks orientation [two-sided unpaired t-test, p = 0.99 (intra-modular half vs. non-module half); p = 0.38 (intra-modular full vs. non-module full)]. In summary, activity decreased as the size of the injury increased, and intra-modular damage had a greater effect than inter-modular damage.
The recoverability of the network was dependent on multiple factors, such as the size and direction of damage, and the presence of modular structure. When the injury was small and applied along the tracks (green bar in Figure 4B), activity recovered to the pre-damage level 2 h later [inter-modular half, 1.04 ± 0.38, two-sided one-sample t-test, p = 0.71 (vs. 1.00)]. This recovery was further observed at 6 h, with activity comparable to that without damage in the cases of damage along tracks (green and red bars) and half damage across tracks (blue bar) (intra-modular half, 0.98 ± 0.57, p = 0.87; inter-modular full, 0.89 ± 0.40, p = 0.32). However, even after 24 h, the activity did not recover in the case of full intra-modular damage (orange bar), i.e., across tracks (intra-modular full, 0.70 ± 0.33, p < 0.05). Furthermore, in the absence of a specific modular structure (purple and pink bars), activity was not restored, regardless of damage size (non-module half, 0.75 ± 0.33, p < 0.05; non-module full, 0.51 ± 0.25, p < 0.001). The differences in damage and recovery caused by network structure and cut orientations were also observed when the same number of connections were eliminated (Supplementary Figure S7). These findings suggest that modular structures facilitate recovery from local injuries.
To examine differences in damage recovery depending on the presence or absence of a modular structure, we focused on community structure before and after damage. We quantified the differences in community structure using normalized mutual information (NMI) (Figure 4C). Immediately after injury, the NMI scores between post- and pre-damage states dropped below 1 for all damage conditions (Figure 4D), indicating that the community structure was altered by the injury. The comparison of the injury sizes revealed that the extent of damage strongly affected the communities in both the intra-modular and non-modular scenarios [two-sided unpaired t-test, p < 0.01 (intra-modular half vs. intra-modular full); p < 0.05 (non-module half vs. non-module full), n = 15]. However, cut size did not affect the communities in the case of inter-modular damage [two-sided unpaired t-test, p = 0.30 (inter-modular half vs. inter-modular full)]. In the context of the tracks pattern, damage across tracks had a more significant impact than damage along the tracks [two-sided unpaired t-test, p < 0.01 (intra-modular half vs. inter-modular half); p < 0.05 (intra-modular half vs. inter-modular full)]. Among these damages on modular networks, the strongest influence was in the networks with full intra-modular damage, but they had greater NMI than networks without a specific modular structure [two-sided unpaired t-test, p < 0.05 (intra-modular full vs. non-module half)]. In summary, we conclude that modular structures reduce the alteration of community structures due to damage, leading to faster recovery than non-modular networks in the SNN.
3.4 STDP restoration of information representation and processing in a spiking neural network model
Finally, we employed a reservoir computing framework to explore the computational impact of network recovery using STDP (Figure 5A). In the framework used here, a signal is injected into the network and its response is linearly decoded to perform a classification task. Spoken digits were used as sensory information and served as input to 5% of the neurons in the SNN both before and after damage. Figure 5B shows the neural responses to the spoken digit “zero” before and after full damage across tracks. Before the injury, the neurons responded well to the input, and the summed activity of the neurons produced two peaks. These peaks of summed neural activity aligned with the peaks of the summed spectrogram of the input spoken digits with a time delay. After injury, the neural response decreased, and the two peaks became smaller. The magnitude of the response recovered after 24 h. This indicates that self-organization through synaptic plasticity restores not only spontaneous activity but also the capacity of the system to respond to external inputs.

Figure 5. Reservoir computing on a damaged spiking neural network. (A) Sketch of the reservoir computing framework. The spectrogram of spoken digits is delivered into a spiking neural network, and the evoked activity is regressed as output. (B) Raster plots of evoked activity for the spoken digit “zero” before (left), immediately after (center), and 24 h after injury (right). The bottom panels show the summed-up spikes (black) and the spectrogram (blue). (C) Time courses of accuracy in the reservoir computing tasks for the six damage conditions. The colors correspond to each type of damage. Gray curves show the undamaged case for the tracks-patterned network, and the brown curves show the undamaged case for the control, unpatterned network. In the panels, the lines indicate the average value, and the shaded areas represent the 95% confidence interval. *p < 0.05; **p < 0.01; ***p < 0.001 (two-sided unpaired t-test, n = 15, df = 13).
Subsequently, we used multiple spoken digits to perform classification tasks. In this task, the linear decoder in the output layer was trained to obtain the optimal output weight matrix Wout to classify the spoken digits “zero,” “one,” and “two.” When the output weight matrix was trained on the network responses before damage, the classification accuracies were 66.2 ± 11.6% (mean ± SD) for the tracks-patterned network and 69.8 ± 9.5% for the unpatterned network, both of which were significantly higher than a random guess (33.3%). When the output weights were fixed at the pre-damage state, classification accuracy decreased over time, regardless of the presence of damage (Supplementary Figure S8). This decline occurred because STDP alters the synaptic weights of the reservoir’s SNN, resulting in changes in the information representation.
Furthermore, we retrained the output weight matrix for each damaged and recovered network state. In the undamaged case, performance was maintained when Wout was retrained (Figure 5C). This indicates that the input information was quantitatively preserved in the neural response and could be classified by a linear decoder, even though the information representation changed over time due to plasticity. Such maintenance of performance was observed under three of the four conditions applied to the modular neural network (intra-modular half, inter-modular half, and inter-modular full), suggesting that the modular structure is robust against small damage in terms of sensory information representation. In contrast, under the conditions of damage on non-module and full intra-modular damage, accuracy dropped immediately after damage (intra-modular full, 59.8 ± 10.1%, p < 0.05; non-module half, 62.2 ± 11.5%, p < 0.05; non-module full, 58.0 ± 11.2%, p < 0.01; n = 15). In other words, under these conditions, damage impaired the network’s ability to represent information and perform speech recognition as a reservoir. While such temporary dysfunction was observed, classification accuracy gradually recovered over time and returned to its original level after 24 h. This suggests that information representation in neural networks improves with the recovery of synchronous activity.
4 Discussion
In the present study, we developed a cultured neuronal network with a modular structure using topographical substrates. Modularity is a distinctive feature of the brain, evolutionarily conserved across many species (Meunier et al., 2010), and is believed to be a crucial trait in living neuronal circuits (Sporns and Betzel, 2016; Michaels et al., 2020; Vishwanathan et al., 2024). Recent studies in neuronal network patterning in vitro have shown that developing a modular structure in neuronal networks can generate rich activity that mimics the complex information processing of the brain (Yamamoto et al., 2018, 2023; Montalà-Flaquer et al., 2022). PDMS-based engineering plays a pivotal role in the development of neuronal cultures with rich structure and dynamics (Halldorsson et al., 2015; Souza et al., 2024). Microfluidic devices made with PDMS can control the structure of neuronal networks (Takemuro et al., 2020), and culturing cortical neurons on PDMS can reduce mechanical mismatch (Sumi et al., 2020) and adjust the topology of the network through topographical modulation (Sharma et al., 2019; Montalà-Flaquer et al., 2022). Furthermore, when damaging neuronal networks, PDMS provides a scaffold that allows mechanical damage to be applied with a scalpel (Ayasreh et al., 2022).
In the present study, we first examined the changes in activity caused by damage to cultured neuronal networks, extending previous studies that also investigated the impact of damage. These studies used laser microdissection or a scalpel to damage a subpopulation of the network (Teller et al., 2020; Ayasreh et al., 2022), while others focused on the modulation of structural connectivity through heat (Hong and Nam, 2020). In these works, the authors observed that cultured neuronal networks temporarily decreased the number of activations and the rate of synchronous activity when damaged, although activity was restored within minutes or days, depending on the magnitude of the damage. However, the overall recovery was not homogeneous; the area of the network surrounding the damage recovered well, whereas the damaged area remained silent and was effectively unrecoverable. Our study showed that immediately after injury, the rate of activity in the network decreased but recovered after 24 h (Figure 2E), which is consistent with previous studies.
Motivated by these culture experiments, we extended our study using a spiking neural network model. The combination of mathematical models and culture experiments facilitated the design of extended in silico explorations, allowing us to compare a wide range of conditions. Specifically, numerical simulations enabled us to explore the mechanisms underlying the restoration of neuronal network activity, which may involve the reorganization of the neuronal network through plasticity. STDP is a natural candidate for such plasticity because it plays important remodeling roles in various regions of the brain (Bi and Poo, 1998; Schmidgall et al., 2024) and has led to the restoration of activity in simulations of randomly inactivating neurons (Gabrieli et al., 2020). Even in cases of severe network disconnection, as in our experiments, STDP worked well and was sufficient to restore neuronal activity (Figure 2J). Further investigation of synaptic weights revealed that synaptic plasticity enhanced global efficiency in the damaged network (Figure 3D). Although previous studies using cultured neuronal networks have estimated the recovery of global efficiency by inferring effective connectivity from network activity (Teller et al., 2020; Ayasreh et al., 2022), our study has shown that it is also possible to examine the synaptic weights underlying firing patterns using a mathematical model that mimics cultured neuronal networks. The qualitative agreement between experiments and simulations opens new avenues for numerically exploring other plasticity mechanisms that are present in cultures, such as homeostasis, and invites experimentalists to devise new experiments to fully monitor synaptic alterations during development and damage.
Related to this, an experimental limitation of our study that needs further investigation is the fact that STDP was not directly demonstrated in the in vitro experiments. This could be achieved for instance through the pharmacological blockade of NMDA receptors in excitatory neurons, whose allosteric kinetics have been related to be central for STDP (Urakubo et al., 2008). While such exploration is promising, it may also affect overall network activity given the NMDA receptors are excitatory, rendering it difficult to isolate the specific contribution of STDP without altering the network dynamics. Experimental design inherently limits the ability to isolate individual contributions within the complex dynamics of living neuronal networks, highlighting the complementary strength of computational models in uncovering the mechanisms that regulate activity in neuronal cultures, specifically in the context of damage and recovery.
The proposed model and its capacity to successfully reproduce the experimental observations allowed us to extend the simulations to explore the effects of different damage conditions. The results showed that the impact of intra-modular damage was stronger than the inter-modular one, and that networks without modular organization were unable to fully recover their activity within the timeframe considered here (Figure 4B). In this regard, we suggest that recovery can be achieved by reorganizing new functional communities that couple regions activating together. Although the new communities departed from the original ones, the underlying initial modular configuration seemed to provide a backbone to facilitate the maintenance of functional communities (Figure 4D). This hypothesis is supported by the fact that modular organization is known to be a damage-resistant structure due to its functional separation and redundancy (Zhang et al., 2019; Chen et al., 2021). For instance, animal experiments have shown that rodents can efficiently perform short-term memory tasks when artificially perturbed by optogenetic stimulation in the brain hemisphere; however, a modular organization is required for the robustness of persistent activity in response to perturbations (Chen et al., 2021). Additionally, in the present study, the burst rate of some numerical explorations, particularly the “half inter-modular” ones (Figure 4B), was significantly higher than in control networks 24 h after injury. Such overactivity is interesting, and may be associated with clinical studies showing that epilepsy can occur after injury (Herman, 2002; Englander et al., 2014; Ding et al., 2016). This overactivity is thought to be triggered by an imbalance between excitation and inhibition. Thus, the wiring of excitatory and inhibitory mesoscale circuits, and their alteration upon damage, may be an aspect that could be explored in a future investigation.
Finally, we observed an association between the restoration of synchronous activity and the recovery of cognitive-like functions, such as spoken digits classification, using a reservoir computing framework. Reservoir computing can be integrated with cultured or computational neural networks to link neuronal dynamics to cognitive and behavioral tasks (Nicola and Clopath, 2017; Yada et al., 2021; Cai et al., 2023; Sumi et al., 2023). In the present study, we used this framework to investigate how changes in neuronal dynamics before and after injury affect speech signal classification. Our model demonstrated that the modular structure was functionally robust to relatively small injuries and that functional impairment, even due to severe damage, could be restored if down-stream neuronal circuits were able to re-organize as well (Figure 5C). Although injury-induced cognitive impairment and recovery have been observed in animal studies (Chen et al., 2010; Christensen et al., 2008; Edlow et al., 2021), our results provide a new perspective on how the dynamic process of synaptic reorganization after injury affects information representation and processing at the neuronal network level. The reservoir computing framework can also be employed for other tasks, such as motor learning and memory, by adjusting the readout. Therefore, it can be applicable to model various neurological dysfunctions and recovery. In our study, functional recovery was achieved through spontaneous synaptic reorganization of functional connectivity and task-dependent learning in the output layer. However, in clinical rehabilitation, patients modify the functional connectivity of the motor cortex through input-output learning in addition to spontaneous modulation (Murata et al., 2015). The recovery through input-output interactions is reported to be promoted by repeated 40 Hz stimulation, which rescues synaptic plasticity such as STDP (Wang et al., 2023). Also, a theoretical study has shown that repeated stimulation inputs to the STDP model can lead to the formation of temporal patterns (Hosaka et al., 2008). The ongoing development of these fields is expected to lead to future applications in rehabilitation models using repetitive input and synaptic plasticity in the neuronal network to restore function.
To conclude, it is worth emphasizing that damage to neuronal networks occurs in several pathologies (Carmichael, 2016; Nagappan et al., 2020). In our study, a cultured neuronal network was designed as an accessible laboratory model for damage to neuronal circuits, which was successfully reproduced computationally. Our model has the potential to predict changes in functional neuronal networks and dynamics due to local damage and may be useful in designing advanced models that combine spontaneous and evoked activity to treat dysfunction. In addition, the mathematical description of the self-organized recovery process can be applied to the field of information processing. The brain is often compared to an electronic computer as an information-processing device. However, damage resistance and recovery capabilities are unique characteristics of living neuronal networks, particularly in the brain (Hassabis et al., 2017). Analysis of recovery from damage through cultured neuronal networks and mathematical neural models is valuable for understanding the mechanisms of neuronal function recovery and for designing artificial neural networks with damage tolerance and self-repair capabilities (Psaier and Dustdar, 2011; Khlaisamniang et al., 2023). Our findings could contribute to further research to introduce self-repair capabilities in robots and AI systems and to advance the understanding of diseases caused by local injury and their treatments in the human nervous system.
Data availability statement
The experimental data are available on Zenodo at https://zenodo.org/records/15314144 (Sumi et al., 2025). The data were analyzed using NETCAL (Orlandi et al., 2017), available at https://github.com/orlandi/netcal. The simulation was based on the Izhikevich polychronization network model (Izhikevich, 2006), available at https://www.izhikevich.org/publications/spnet.htm. The transfer entropy analysis was performed using code from https://github.com/orlandi/te-causality, and the axonal growth model used to generate the network structure is available at https://github.com/akkeh/nc_sim.
Ethics statement
The animal study was approved by the Ethics Committee for Animal Experimentation of the University of Barcelona (CEEA-UB). The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
TS: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Writing – original draft, Writing – review & editing. AH: Formal Analysis, Methodology, Validation, Writing – review & editing, Software. HY: Conceptualization, Funding acquisition, Supervision, Validation, Writing – original draft, Writing – review & editing. HK: Methodology, Validation, Writing – review & editing, Software. YK: Writing – review & editing, Funding acquisition, Conceptualization, Validation. JS: Conceptualization, Data curation, Funding acquisition, Methodology, Supervision, Validation, Writing – original draft, Writing – review & editing, Software. AH-I: Supervision, Validation, Writing – review & editing, Funding acquisition.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by the MEXT Grant-in-Aid for Transformative Research Areas (A) and (B) “Multicellular Neurobiocomputing” (21H05163, 21H05164, 24H02330, 24H02332, 24H02334), JSPS KAKENHI (21K12050, 22KK0177, 23H00251, 22KJ0190, 23K28179), the JST Moonshot R&D Program (JPMJMS2023-41), the JSPS Overseas Challenge Program for Young Researchers, and the Tohoku University Research Institute of Electrical Communication (RIEC) Cooperative Research Project Program. This research was supported by the European Union Horizon 2020 research and innovation program under grant 964877 (project NEU-CHiP). JS acknowledged support from grant PID2022-137713NB-C22, funded by MCIU/AEI/10.13039/501100011033, ERDF/EU, and the Generalitat de Catalunya under grant 2021-SGR-00450.
Acknowledgments
We are grateful to Hayato Chiba for valuable discussions.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The authors declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2025.1570783/full#supplementary-material
References
Alotaibi, N., and Rhouma, D. (2022). A review on community structures detection in time evolving social networks. J. King Saud Univ. Comput. Inform. Sci. 34, 5646–5662. doi: 10.1016/j.jksuci.2021.08.016
Arnemann, K. L., Chen, A. J.-W. Novakovic-Agopian, T., Gratton, C., and Nomura, E. M.D’Esposito, M. (2015). Functional brain network modularity predicts response to cognitive training after brain injury. Neurology 84, 1568–1574. doi: 10.1212/WNL.0000000000001476
Ayasreh, S., Jurado, I., López-León, C. F., Montalà-Flaquer, M., and Soriano, J. (2022). Dynamic and functional alterations of neuronal networks in vitro upon physical damage: A proof of concept. Micromachines 13:2259. doi: 10.3390/mi13122259
Bi, G., and Poo, M. (1998). Synaptic modifications in cultured hippocampal neurons: Dependence on spike timing, synaptic strength, and postsynaptic cell type. J. Neurosci. 18, 10464–10472. doi: 10.1523/JNEUROSCI.18-24-10464.1998
Blondel, V. D., Guillaume, J.-L., Lambiotte, R., and Lefebvre, E. (2008). Fast unfolding of communities in large networks. J. Stat. Mech. 2008:10008. doi: 10.1088/1742-5468/2008/10/P10008
Boroda, E., Armstrong, M., Gilmore, C. S., Gentz, C., Fenske, A., Fiecas, M., et al. (2021). Network topology changes in chronic mild traumatic brain injury (mTBI). Neuro. Image Clinical. 31:102691. doi: 10.1016/j.nicl.2021.102691
Bunday, K. L., and Perez, M. A. (2012). Motor recovery after spinal cord injury enhanced by strengthening corticospinal synaptic transmission. Curr. Biol. 22, 2355–2361. doi: 10.1016/j.cub.2012.10.046
Cai, H., Ao, Z., Tian, C., Wu, Z., Liu, H., Tchieu, J., et al. (2023). Brain organoid reservoir computing for artificial intelligence. Nat. Electron. 6, 1032–1039. doi: 10.1038/s41928-023-01069-w
Carmichael, S. T. (2016). The 3 Rs of stroke biology: Radial, relayed, and regenerative. Neurotherapeutics 13, 348–359. doi: 10.1007/s13311-015-0408-0
Chen, G., Kang, B., Lindsey, J., Druckmann, S., and Li, N. (2021). Modularity and robustness of frontal cortical networks. Cell 184, 3717–3730.e24. doi: 10.1016/j.cell.2021.05.026.
Chen, H., Epstein, J., and Stern, E. (2010). Neural plasticity after acquired brain injury: Evidence from functional neuroimaging. PM&R 2, S306–S312. doi: 10.1016/j.pmrj.2010.10.006
Christensen, B. K., Colella, B., Inness, E., Hebert, D., Monette, G., Bayley, M., et al. (2008). Recovery of cognitive function after traumatic brain injury: A multilevel modeling analysis of Canadian outcomes. Arch. Phys. Med. Rehabil. 89, S3–S15. doi: 10.1016/j.apmr.2008.10.002
Ding, K., Gupta, P. K., and Diaz-Arrastia, R. (2016). “Epilepsy after traumatic brain injury,” in Translational Research in Traumatic Brain Injury, eds D. Laskowitz and G. Grant (Boca Raton, FL: CRC Press).
Edlow, B. L., Claassen, J., Schiff, N. D., and Greer, D. M. (2021). Recovery from disorders of consciousness: Mechanisms, prognosis and emerging therapies. Nat. Rev. Neurol. 17, 135–156. doi: 10.1038/s41582-020-00428-x
Effenberger, F., Jost, J., and Levina, A. (2015). Self-organization in balanced state networks by STDP and homeostatic plasticity. PLoS Comput. Biol. 11:e1004420. doi: 10.1371/journal.pcbi.1004420
Englander, J., Cifu, D. X., and Diaz-Arrastia, R. (2014). Seizures after traumatic brain injury. Arch. Phys. Med. Rehabil. 95, 1223–1224. doi: 10.1016/j.apmr.2013.06.002
Fagiolo, G. (2007). Clustering in complex directed networks. Phys. Rev. E 76:026107. doi: 10.1103/PhysRevE.76.026107
Fukui, A., Osaki, H., Ueta, Y., Kobayashi, K., Muragaki, Y., Kawamata, T., et al. (2020). Layer-specific sensory processing impairment in the primary somatosensory cortex after motor cortex infarction. Sci. Rep. 10:3771. doi: 10.1038/s41598-020-60662-7
Gabrieli, D., Schumm, S. N., Vigilante, N. F., Parvesse, B., and Meaney, D. F. (2020). Neurodegeneration exposes firing rate dependent effects on oscillation dynamics in computational neural networks. PLoS One 15:e0234749. doi: 10.1371/journal.pone.0234749
Gütig, R., Aharonov, R., Rotter, S., and Sompolinsky, H. (2003). Learning input correlations through nonlinear temporally asymmetric hebbian plasticity. J. Neurosci. 23, 3697–3714. doi: 10.1523/JNEUROSCI.23-09-03697.2003
Halldorsson, S., Lucumi, E., Gómez-Sjöberg, R., and Fleming, R. M. T. (2015). Advantages and challenges of microfluidic cell culture in polydimethylsiloxane devices. Biosens. Bioelectron. 63, 218–231. doi: 10.1016/j.bios.2014.07.029
Hassabis, D., Kumaran, D., Summerfield, C., and Botvinick, M. (2017). Neuroscience-inspired artificial intelligence. Neuron 95, 245–258. doi: 10.1016/j.neuron.2017.06.011
Herman, S. T. (2002). Epilepsy after brain insult. Neurology 59, S21–S26. doi: 10.1212/WNL.59.9_suppl_5.S21
Hong, N., and Nam, Y. (2020). Thermoplasmonic neural chip platform for in situ manipulation of neuronal connections in vitro. Nat. Commun. 11:6313. doi: 10.1038/s41467-020-20060-z
Hosaka, R., Araki, O., and Ikeguchi, T. (2008). STDP provides the substrate for igniting synfire chains by spatiotemporal input patterns. Neural Comput. 20, 415–435. doi: 10.1162/neco.2007.11-05-043
Houben, A. M., Garcia-Ojalvo, J., and Soriano, J. (2025). Role of connectivity anisotropies in the dynamics of cultured neuronal networks. arXiV [Preprint] doi: 10.48550/arXiv.2501.04427
Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14, 1569–1572. doi: 10.1109/TNN.2003.820440
Izhikevich, E. M. (2006). Polychronization: Computation with spikes. Neural Comput. 18, 245–282. doi: 10.1162/089976606775093882
Izhikevich, E. M., and Desai, N. S. (2003). Relating STDP to BCM. Neural Comput. 15, 1511–1523. doi: 10.1162/089976603321891783
Khlaisamniang, P., Khomduean, P., Saetan, K., and Wonglapsuwan, S. (2023). “Generative AI for self-healing systems,” in Proceedings of the 2023 18th International Joint Symposium on Artificial Intelligence and Natural Language Processing (iSAI-NLP), (Bangkok: Curran Associates, Inc), 1–6. doi: 10.1109/iSAI-NLP60301.2023.10354608
Latora, V., and Marchiori, M. (2001). Efficient behavior of small-world networks. Phys. Rev. Lett. 87:198701. doi: 10.1103/PhysRevLett.87.198701
Lee, W.-C. A., Bonin, V., Reed, M., Graham, B. J., Hood, G., Glattfelder, K., et al. (2016). Anatomy and function of an excitatory network in the visual cortex. Nature 532, 370–374. doi: 10.1038/nature17192
Liberman, M., Amsler, R., Church, K., Fox, E., Hafner, C., Klavans, J., et al. (1993). TI 46-Word. Philadelphia: Linguistic Data Consortium, doi: 10.35111/ZX7A-FW03
Lynn, C. W., and Bassett, D. S. (2019). The physics of brain network structure, function and control. Nat. Rev. Phys. 1, 318–332. doi: 10.1038/s42254-019-0040-8
Meunier, D., Lambiotte, R., and Bullmore, E. T. (2010). Modular and hierarchically modular organization of brain networks. Front. Neurosci. 4:200. doi: 10.3389/fnins.2010.00200
Michaels, J. A., Schaffelhofer, S., Agudelo-Toro, A., and Scherberger, H. (2020). A goal-driven modular neural network predicts parietofrontal neural dynamics during grasping. Proc. Natl. Acad. Sci. U. S. A. 117, 32124–32135. doi: 10.1073/pnas.2005087117
Montalà-Flaquer, M., López-León, C. F., Tornero, D., Houben, A. M., Fardet, T., Monceau, P., et al. (2022). Rich dynamics and functional organization on topographically designed neuronal networks in vitro. iScience 25:105680. doi: 10.1016/j.isci.2022.105680
Morrison, A., Diesmann, M., and Gerstner, W. (2008). Phenomenological models of synaptic plasticity based on spike timing. Biol. Cybern 98, 459–478. doi: 10.1007/s00422-008-0233-1
Murata, Y., Higo, N., Hayashi, T., Nishimura, Y., Sugiyama, Y., Oishi, T., et al. (2015). Temporal plasticity involved in recovery from manual dexterity deficit after motor cortex lesion in macaque monkeys. J. Neurosci. 35, 84–95. doi: 10.1523/JNEUROSCI.1737-14.2015
Nagappan, P. G., Chen, H., and Wang, D.-Y. (2020). Neuroregeneration and plasticity: A review of the physiological mechanisms for achieving functional recovery postinjury. Military Med. Res. 7:30. doi: 10.1186/s40779-020-00259-3
Nicola, W., and Clopath, C. (2017). Supervised learning in spiking neural networks with FORCE training. Nat. Commun. 8:2208. doi: 10.1038/s41467-017-01827-3
Nudo, R. J. (2013). Recovery after brain injury: Mechanisms and principles. Front. Hum. Neurosci. 7:887. doi: 10.3389/fnhum.2013.00887
O’Neill, K. M., Anderson, E. D., Mukherjee, S., Gandu, S., McEwan, S. A., Omelchenko, A., et al. (2023). Time-dependent homeostatic mechanisms underlie brain-derived neurotrophic factor action on neural circuitry. Commun. Biol. 6, 1–26. doi: 10.1038/s42003-023-05638-9
Onnela, J.-P., Saramäki, J., Kertész, J., and Kaski, K. (2005). Intensity and coherence of motifs in weighted complex networks. Phys. Rev. E 71:065103. doi: 10.1103/PhysRevE.71.065103
Orlandi, J. G., Fernández-García, S., Comella-Bolla, A., Masana, M., Barriga, G. G.-D., Yaghoubi, M., et al. (2017). NETCAL: An Interactive Platform for Large-Scale, NETwork and Population Dynamics Analysis of Calcium Imaging Recordings. Switzerland: Zenoda. doi: 10.5281/zenodo.1119026
Orlandi, J. G., Soriano, J., Alvarez-Lacalle, E., Teller, S., and Casademunt, J. (2013). Noise focusing and the emergence of coherent activity in neuronal cultures. Nat. Phys. 9, 582–590. doi: 10.1038/nphys2686
Pasquale, V., Martinoia, S., and Chiappalone, M. (2017). Stimulation triggers endogenous activity patterns in cultured cortical networks. Sci. Rep. 7:9080. doi: 10.1038/s41598-017-08369-0
Psaier, H., and Dustdar, S. (2011). A survey on self-healing systems: Approaches and systems. Computing 91, 43–73. doi: 10.1007/s00607-010-0107-y
Rubin, J., Lee, D. D., and Sompolinsky, H. (2001). Equilibrium properties of temporally asymmetric hebbian plasticity. Phys. Rev. Lett. 86, 364–367. doi: 10.1103/PhysRevLett.86.364
Rubinov, M., and Sporns, O. (2010). Complex network measures of brain connectivity: Uses and interpretations. Neuroimage 52, 1059–1069. doi: 10.1016/j.neuroimage.2009.10.003
Schmidgall, S., Ziaei, R., Achterberg, J., Kirsch, L., Hajiseyedrazi, S. P., and Eshraghian, J. (2024). Brain-inspired learning in artificial neural networks: A review. APL Mach. Learn. 2:021501. doi: 10.1063/5.0186054
Sharma, D., Jia, W., Long, F., Pati, S., Chen, Q., Qyang, Y., et al. (2019). Polydopamine and collagen coated micro-grated polydimethylsiloxane for human mesenchymal stem cell culture. Bioactive Mater. 4, 142–150. doi: 10.1016/j.bioactmat.2019.02.002
Slaney, M., and Lyon, R. (1993). On the Importance of time—a Temporal Representation of Sound. Available online at: https://www.semanticscholar.org/paper/On-the-importance-of-time%E2%80%94a-temporal-representation-Slaney-Lyon/cbfc0cc32d1670a2c5421b2872fa6afadcb796b9 (accessed November 26, 2024).
Soriano, J. (2023). Neuronal cultures: Exploring biophysics, complex systems, and medicine in a dish. Biophysica 3, 181–202. doi: 10.3390/biophysica3010012
Souza, A., Nobrega, G., Neves, L. B., Barbosa, F., Ribeiro, J., Ferrera, C., et al. (2024). Recent Advances of PDMS in vitro biomodels for flow visualizations and measurements: From macro to nanoscale applications. Micromachines 15:1317. doi: 10.3390/mi15111317
Sporns, O., and Betzel, R. F. (2016). Modular brain networks. Annu. Rev. Psychol. 67, 613–640. doi: 10.1146/annurev-psych-122414-033634
Stetter, O., Battaglia, D., Soriano, J., and Geisel, T. (2012). Model-free reconstruction of excitatory neuronal connectivity from calcium imaging signals. PLoS Comput. Biol. 8:e1002653. doi: 10.1371/journal.pcbi.1002653
Sumi, T., Houben, A. M., Yamamoto, H., Kato, H., Katori, Y., Soriano Fradera, J., et al. (2025). Dataset for: Modular architecture confers robustness to damage and facilitates recovery in spiking neural networks modeling in vitro neurons. Front. Neurosci. 19:15314144. doi: 10.5281/zenodo.15314144
Sumi, T., Yamamoto, H., and Hirano-Iwata, A. (2020). Suppression of hypersynchronous network activity in cultured cortical neurons using an ultrasoft silicone scaffold. Soft. Matter. 16, 3195–3202. doi: 10.1039/C9SM02432H
Sumi, T., Yamamoto, H., Katori, Y., Ito, K., Moriya, S., Konno, T., et al. (2023). Biological neurons act as generalization filters in reservoir computing. Proc. Natl. Acad. Sci. 120:e2217008120. doi: 10.1073/pnas.2217008120
Takemuro, T., Yamamoto, H., Sato, S., and Hirano-Iwata, A. (2020). Polydimethylsiloxane microfluidic films for in vitro engineering of small-scale neuronal networks. Jpn. J. Appl. Phys. 59:117001. doi: 10.35848/1347-4065/abc1ac
Teller, S., Estévez-Priego, E., Granell, C., Tornero, D., Andilla, J., Olarte, O. E., et al. (2020). Spontaneous functional recovery after focal damage in neuronal cultures. eNeuro 7:ENEURO.0254-19.2019. doi: 10.1523/ENEURO.0254-19.2019
Tibau, E., Ludl, A.-A., Rüdiger, S., Orlandi, J. G., and Soriano, J. (2020). Neuronal spatial arrangement shapes effective connectivity traits of in vitro cortical networks. IEEE Trans. Netw. Sci. Eng. 7, 435–448. doi: 10.1109/TNSE.2018.2862919
Turrigiano, G. G. (1999). Homeostatic plasticity in neuronal networks: The more things change, the more they stay the same. Trends Neurosci. 22, 221–227. doi: 10.1016/S0166-2236(98)01341-1
Turrigiano, G. G. (2008). The self-tuning neuron: Synaptic scaling of excitatory synapses. Cell 135, 422–435. doi: 10.1016/j.cell.2008.10.008
Urakubo, H., Honda, M., Froemke, R. C., and Kuroda, S. (2008). Requirement of an allosteric kinetics of NMDA receptors for spike timing-dependent plasticity. J. Neurosci. 28, 3310–3323. doi: 10.1523/JNEUROSCI.0303-08.2008
van Rossum, M. C. W., Bi, G. Q., and Turrigiano, G. G. (2000). Stable hebbian learning from spike timing-dependent plasticity. J. Neurosci. 20, 8812–8821. doi: 10.1523/JNEUROSCI.20-23-08812.2000
Vishwanathan, A., Sood, A., Wu, J., Ramirez, A. D., Yang, R., Kemnitz, N., et al. (2024). Predicting modular functions and neural coding of behavior from a synaptic wiring diagram. Nat. Neurosci. 27, 2443–2454. doi: 10.1038/s41593-024-01784-3
Wagenaar, D. A., Pine, J., and Potter, S. M. (2006). An extremely rich repertoire of bursting patterns during the development of cortical cultures. BMC Neurosci. 7:11. doi: 10.1186/1471-2202-7-11
Walker, K. R., and Tesco, G. (2013). Molecular mechanisms of cognitive dysfunction following traumatic brain injury. Front. Aging Neurosci. 5:29. doi: 10.3389/fnagi.2013.00029
Wang, C., Lin, C., Zhao, Y., Samantzis, M., Sedlak, P., Sah, P., et al. (2023). 40-Hz optogenetic stimulation rescues functional synaptic plasticity after stroke. Cell Rep. 42:113475. doi: 10.1016/j.celrep.2023.113475
Yada, Y., Yasuda, S., and Takahashi, H. (2021). Physical reservoir computing with FORCE learning in a living neuronal culture. Appl. Phys. Lett. 119:173701. doi: 10.1063/5.0064771
Yamamoto, H., Kubota, S., Chida, Y., Morita, M., Moriya, S., Akima, H., et al. (2016). Size-dependent regulation of synchronized activity in living neuronal networks. Phys. Rev. E 94:012407. doi: 10.1103/PhysRevE.94.012407
Yamamoto, H., Moriya, S., Ide, K., Hayakawa, T., Akima, H., Sato, S., et al. (2018). Impact of modular organization on dynamical richness in cortical networks. Sci. Adv. 4:eaau4914. doi: 10.1126/sciadv.aau4914
Yamamoto, H., Spitzner, F. P., Takemuro, T., Buendía, V., Murota, H., Morante, C., et al. (2023). Modular architecture facilitates noise-driven control of synchrony in neuronal networks. Sci. Adv. 9:eade1755. doi: 10.1126/sciadv.ade1755
Keywords: spiking neural network, spike timing dependent plasticity, cultured neuronal network, modular structure, reservoir computing
Citation: Sumi T, Houben AM, Yamamoto H, Kato H, Katori Y, Soriano J and Hirano-Iwata A (2025) Modular architecture confers robustness to damage and facilitates recovery in spiking neural networks modeling in vitro neurons. Front. Neurosci. 19:1570783. doi: 10.3389/fnins.2025.1570783
Received: 14 February 2025; Accepted: 20 May 2025;
Published: 19 June 2025.
Edited by:
Michela Chiappalone, University of Genoa, ItalyReviewed by:
Paolo Massobrio, University of Genoa, ItalyFikret Emre Kapucu, Tampere University, Finland
Copyright © 2025 Sumi, Houben, Yamamoto, Kato, Katori, Soriano and Hirano-Iwata. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Takuma Sumi, dGFrdW1hLnN1bWkuYzhAdG9ob2t1LmFjLmpw