Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Reaction-diffusion models in weighted and directed connectomes

  • Oliver Schmitt ,

    Contributed equally to this work with: Oliver Schmitt, Christian Nitzsche

    Roles Conceptualization, Data curation, Formal analysis, Investigation, Resources, Supervision, Visualization, Writing – original draft, Writing – review & editing

    oliver.schmitt@medicalschool-hamburg.de, schmitt@med.uni-rostock.de

    Affiliation Medical School Hamburg, University of Applied Sciences, Hamburg; University of Rostock, Department of Anatomy, Rostock, Germany

    ⨯
  • Christian Nitzsche ,

    Contributed equally to this work with: Oliver Schmitt, Christian Nitzsche

    Roles Methodology, Software, Validation

    Affiliation Medical School Hamburg, University of Applied Sciences, Hamburg; University of Rostock, Department of Anatomy, Rostock, Germany

    ⨯
  • Peter Eipert,

    Roles Software, Validation

    Affiliation Medical School Hamburg, University of Applied Sciences, Hamburg; University of Rostock, Department of Anatomy, Rostock, Germany

    ⨯
  • Vishnu Prathapan,

    Roles Software, Validation, Writing – original draft

    Affiliation Medical School Hamburg, University of Applied Sciences, Hamburg; University of Rostock, Department of Anatomy, Rostock, Germany

    ⨯
  • Marc-Thorsten Hütt,

    Roles Conceptualization, Methodology, Writing – original draft

    Affiliation Jacobs University, Life Sciences & Chemistry, Computational Systems Biology, Bremen, Germany

    ⨯
  • Claus C. Hilgetag

    Roles Conceptualization, Writing – original draft

    Affiliation University of Hamburg, Department of Computational Neuroscience, Hamburg, Germany

    ⨯

Abstract

Connectomes represent comprehensive descriptions of neural connections in a nervous system to better understand and model central brain function and peripheral processing of afferent and efferent neural signals. Connectomes can be considered as a distinctive and necessary structural component alongside glial, vascular, neurochemical, and metabolic networks of the nervous systems of higher organisms that are required for the control of body functions and interaction with the environment. They are carriers of functional phenomena such as planning behavior and cognition, which are based on the processing of highly dynamic neural signaling patterns. In this study, we examine more detailed connectomes with edge weighting and orientation properties, in which reciprocal neuronal connections are also considered. Diffusion processes are a further necessary condition for generating dynamic bioelectric patterns in connectomes. Based on our precise connectome data, we investigate different diffusion-reaction models to study the propagation of dynamic concentration patterns in control and lesioned connectomes. Therefore, differential equations for modeling diffusion were combined with well-known reaction terms to allow the use of connection weights, connectivity orientation and spatial distances.

Three reaction-diffusion systems Gray-Scott, Gierer-Meinhardt and Mimura-Murray were investigated. For this purpose, implicit solvers were implemented in a numerically stable reaction-diffusion system within the framework of neuroVIISAS. The implemented reaction-diffusion systems were applied to a subconnectome which shapes the mechanosensitive pathway that is strongly affected in the multiple sclerosis demyelination disease. It was found that demyelination modeling by connectivity weight modulation changes the oscillations of the target region, i.e. the primary somatosensory cortex, of the mechanosensitive pathway.

In conclusion, a new application of reaction-diffusion systems to weighted and directed connectomes has been realized. Because the implementation was realized in the neuroVIISAS framework many possibilities for the study of dynamic reaction-diffusion processes in empirical connectomes as well as specific randomized network models are available now.

Author summary

Reaction-diffusion systems were adapted and analyzed in weighted and directed connectomes. The systems were applied to a multiple sclerosis model by modulating connectivity weights within the reaction-diffusion process. This leads to changes in the oscillation patterns of a target region of the mechanosensitive pathway.

Citation: Schmitt O, Nitzsche C, Eipert P, Prathapan V, Hütt M-T, Hilgetag CC (2022) Reaction-diffusion models in weighted and directed connectomes. PLoS Comput Biol 18(10): e1010507. https://doi.org/10.1371/journal.pcbi.1010507

Editor: Marcus Kaiser, University of Nottingham, UNITED KINGDOM

Received: March 3, 2022; Accepted: August 22, 2022; Published: October 28, 2022

Copyright: © 2022 Schmitt et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the manuscript and its Supporting information files. The source code is available in the public repository figshare and can be downloaded from the web address: https://doi.org/10.6084/m9.figshare.21081418.v1.

Funding: OS received funding by Novartis (MOMB157A_FVAA004). https://www.novartis.com/ The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Diffusion is the process by which matter or particles like atoms or molecules naturally move from regions where they are highly concentrated to regions where they are not as concentrated. We can think of diffusion of different particles as a dynamic process along a line or one dimension. The particles move randomly, which is also called a random walk. In view of the numerous realizations of diffusion processes in natural biological systems and compartments, it appears obvious to apply diffusion in neural networks. Models of reaction-diffusion (RD) adapted to networks can be used to investigate the spreading of information through networks or connectomes which may shape dynamic states that can be related to functional processes. It has been shown that dynamic models of brain communication have begun to create links between connectional architectures and function. Furthermore, brains have the capacity to support a great diversity of dynamic patterns which are complex at a broad range of temporal frames to sustain a large number of competing functional demands. A large amount of different dynamic patterns has been considered as a functional repertoire of a network that allows flexibility across a broad range of cognitive functions [1].

Applications of diffusion

Reaction-diffusion systems produce complex self-organized patterns, such as spreading pulses and fronts, stationary dissipative shapes, rotating waves and turbulences [2–6]. Some reaction-diffusion models have been found to work on networks through interacting species occupying network nodes and diffusively transferred across links [7–11]. Reaction-diffusion models can also correspond to networks of diffusively coupled biological cells or regions and chemical reactors [8–11].

Reaction-diffusion processes can be considered and applied on different scales. On the macro and meso level, the application is suitable for neuronal networks, connectomes and local circuits. Recently, a paper on 3D reaction-diffusion models of neurons and networks was published [12]. At the micro level, reaction-diffusion systems are applied to the modeling of cellular and subcellular processes. Here, NeuroRD (https://github.com/neurord; http://modeldb.science/modellist?id=139757&all_simu=true) should be mentioned, which is used for the simulation of neuronal signaling pathways such as that of stochastic simulation on spiny dendrites [13]. In NeuroRD the Gillespie’s tau-leap reaction algorithm, and the stochastic diffusion algorithm of Blackwell are realized. At the synaptic level [14] and on biological surfaces [15] diffusion modeling is applied, for example, by MCell (https://mcell.org/). In detailed synpase models, diffusion processes are also considered [16, 17]. Beyond MCell, diffusion processes are also used in particularly detailed synpase models [16, 17]. Virtual Cell (VCell) (https://vcell.org/) is a comprehensive platform for modeling cell biological systems based on a central database and distributed as a web application [18, 19]. This environment is also used to model cellular diffusion processes [20, 21] and neuronal signaling [22]. A set of libraries for GENESIS is Chemesis, which allows the simulation of reaction-diffusion systems, including calcium release. This approach has been used by [23] to create in a multicompartment-model of phototransduction, calcium dynamics, and ionic currents of a photoreceptor. We should also mention that there is a wide range of applications for electrodiffusion such as the calculation of extracellular electric potentials from neuron stimulations [24], the simulation of electrodiffusion and water movements in brain tissue [25], ionic electrodiffusion with cortical propagation depression [26], and application in electrodiffusive neuron-extracellular-glia models to study the emergence of slow potentials in the brain [27]. Finally, we must mention here the application of diffusion processes to the modeling of diseases such as strokes [28, 29] and natural developmental processes [30].

In complex networks, the analysis of self-organization is difficult and analysis has been restricted to non-equilibrium pattern formation as synchronization [31–33] or epidemic spreading [34–36]. Turing demonstrated [2] that changes of the diffusion constants of activators and inhibitors produce a destabilization of the uniform state of a system and lead to the spontaneous emergence of oscillatory patterns (Turing patterns) which emerge in chemical reactions, biological morphogenesis, and in ecosystems. This kind of complex non-equilibrium self-organization and Turing instabilities can occur in networks as well [8]. Such network-instabilities have been explored further in mathematical frameworks [9–11]. However, their analyses were restricted to regular lattices [8, 9] and small networks [10, 11]. Recent theoretical work has elucidated the relationship between network architecture and diversity of Turing patterns [37].

Diffusion models have been proven useful for delineating functional modules [38, 39] and predicting statistical dependencies (functional connectivity) among remote neuronal time courses [40, 41]. Network diffusion or graph diffusion [42] has also been applied successfully to network modularity analysis [43] as well as to modeling the relation between structural and functional brain connectivity networks [44].

Meanwhile, several research groups are engaged in the application of diffusion models in networks of nervous systems [40, 44–46]. In most cases, the above-mentioned particles of diffusion processes are considered as components of information. Diffusion processes can then be extended as a flow-based communication model between areas of a nervous system [47–51]. This type of dynamic model is also suitable to describe the phenomenon of multisensory integration [47] or to model the dynamics of brain diseases [45, 46, 48].

Diffusion in connectomes

Through the use of different reaction terms in the differential equations, reaction-diffusion models allow a versatile parameterization of the dynamic behavior in coupled systems such as neuronal connectomes. In the present study, we adopted different models such as predator-prey and activator-inhibitor for use in connectomes to have more options available for studying the propagation of diffusion processes in networks.

Analyzing reaction-diffusion models on connectome architectures provides insight in the pattern forming capabilities and, hence, the feasible collective modes, of such architectures. Here we first illustrate, using simple, generic network architectures, how reaction-diffusion systems create sets of nodes with common dynamical behaviors, which cannot be trivially derived from the network architecture alone. We subsequently apply this approach to the spinal cord, brainstem, diencephalic and cortical connectivity of the mechanosensory pathways. This new approach—probing connectomes with reaction-diffusion models—is fully integrated in neuroVIISAS [52]. A detailed tutorial is provided as Supporting Information.

For the application of reaction-diffusion models, the weights of neuronal connections of a connectome are considered as strengths of connections between regions. Weights of connections in tract-tracing studies encode the number of nerve fibers that have incorporated a tract-tracing substance or an estimate of the number of traced nerve fibers. Thus, estimates of weights of connections do not describe mean thickness of myelin sheaths which are investigated by transmission electron microscopy or other techniques. In most tract-tracing studies weights or densities of connections are estimates which are described in the three basic categories weak, moderate and strong. In addition, connections may be described without any categories at all or in further categories of weights like “weak to moderate” or “very strong”. A survey of all categories, their interpretation, relations and comparisons are given by Schwanke et al. [53]. In terms of a relation of weights of connections and their functional mean they can be interpreted as the strength of a connection which allows transmission of an electrochemical signal between regions of the nervous system. If the strength of a connection is strong then more information can be transmitted. With regard to neuropathological changes of nerve fibers like those in multiple sclerosis where demyelination and axon degeneration occur, the weights of connection can be related to such pathological processes. In this case the weights of connection will be reduced. If weights of connections are used in this study, the ordinal scaled estimates were always logarithmically transformed.

A generalized form of the predator-prey model of Lotka-Volterra was introduced by Mimura-Murray [6, 54, 55] with some advantages with regard to the original concept. The Mimura-Murray model (MM) appears to be promising for adaption to network diffusion.

An activator-inhibitor reaction-diffusion model was developed by Gierer and Meinhardt (GM) [56–58]. It consists of a reaction term with activator and inhibitor parameters that could be adopted to generate specific oscillation patterns in weighted digraphs. Therefore, the GM model appears to be an interesting candidate for applying in network diffusion. Finally, the Gray-Scott model (GS), which is a classical mathematical model for isothermal autocatalytic reaction with another type of reaction term in the differential equation [59–65], has been adapted to network diffusion.

The objective of this study was to investigate reaction-diffusion models with regard to weighted digraphs and distances in the diffusion terms. Here, data of the connectome of the rat nervous system were used and synthetic randomized directed networks with preserved edges and nodes. The reaction-diffusion models of Gierer-Meinhardt [56–58], Mimura-Murray [6, 54, 55] and Gray-Scott [59–65] were adapted to weighted digraphs under consideration of the estimated or linear Euclidean distances between nodes in order to build more realistic coupled dynamic models based on empirical data.

A second aim of this investigation was the analysis of pattern forming in such a way that the reaction-diffusion models generate sets of nodes with common dynamical behavior.

A further objective of this work was to investigate the effects of changes of connection weights within the reaction-diffusion process because this could be a starting point for modeling progression of neurodegenerative diseases [66–71]. The change in connection weights in parameterized models relates specifically to the class of neurodegenerative diseases found in neuronal connections rather than gray matter or neuronal perikarya. Multiple sclerosis is one of such demyelinating diseases [72–81] with a large variety of temporal and topographic progression patterns. The framework used here allows to define sets of source nodes or regions and sets of target nodes of interest embedded in a large connectome context to systematically investigate signal propagation or information diffusion through highly complex connectional architectures which undergo precisely defined changes of connections weights. Hence, this simulation environment appears to be a promising starting point to investigate in a consistent and reproducible way changes of network features as in neurological disorders and their effect of dynamic pattern modifications.

Materials and methods

Connectome data and structure

The connectome data were generated in a metastudy of original research publications. Collating information of neuronal connections between pairs of regions can be performed by manual readout of data in over 7000 original research publications which describe the anterograde and retrograde transport of tract-tracing substances. This metastudy approach is well established and has been successfully performed in ferret, avian, macaque, cat and rat [82–87]. The connectome, circuit and lattice data used in this investigation can be downloaded from [88] and the neuroVIISAS framework from https://neuroviisas.med.uni-rostock.de/neuroviisas.shtml. How to install and start neuroVIISAS is described in the supplement. The test.brain and MS.brain project files can be directly loaded in neuroVIISAS. The reference.bib file (included in [88]) provides the links of connections with original research literature and need to be configured in neuroVIISAS [52].

Implementation of reaction-diffusion models

Starting from the broad application of reaction-diffusion systems presented in the introduction, we would like to investigate the property of pulse propagations in networks. More specifically, the reaction-diffusion models will be applied to weighted and directed connectomes. Since all areas are represented in our connectome data that are also affected in multiple sclerosis diseases, a differential investigation of the reaction-diffusion models of control connectome and in multiple sclerosis lesioned connectome will be performed. Since there are no studies on the propagation of concentrations of the reaction-diffusion models in a directional and weighted connectome, we have implemented three basic models, namely the Gierer-Meinhardt, Mimura-Murray and Gray-Scott models, in the neuroVIISAS framework in such a way that the models are directly applicable with selectable parameters to connectomes with the properties listed above.

In the tutorial part 2 the Gierer-Meinhardt model, as well as the other two models, will be explained in more detail, therefore only the functions of the two substances U(x, y, t) and V(x, y, t) at node N(x, y) of two-dimensional regular lattice should be shown in the following: (1) (2)

The Laplace operator Δ in the following form (3) allows the diffusion of the two substances to be determined.

The reaction-diffusion models applied to the connectome data are based on differential equation systems which could be stiff and are solved by implicit solvers. Their form, solution and extension with regard to noise, connection weight modulation and Euclidean distances are described in detail in the second part of the Supporting Information (Tutorial part 2).

In tutorial part 2, we describe in detail the motivation for using reaction-diffusion models in directed and weighted connectomes. In particular, the formal implementation of the reaction-diffusion models is discussed in detail. For this purpose we start from the Laplace operator to model spatial diffusion. First, the model is not developed for connectomes but for regular lattices. By means of discretization we obtain a system of ordinary differential equations (ODE). For the solution of the ODEs different methods have been implemented in neuroVIISAS (Euler solver, Runge-Kutta, step-controlled Dormand-Prince).

The three reaction-diffusion models are formally described and the initial values of the parameters of the models with which the diffusion patterns known from the literature can be generated. With our implemented models, the same patterns as the patterns known from the literature have been computed, so that our models provide reproducible results. After the formal introduction of the three models, we describe in the tutorial how they are incorporated in a directed network and what has to be considered to define the diffusion directions in a directed network. For this purpose we have provided a selection option that allows to choose between input diffusion and output diffusion of the adjacency matrix and the adjacency matrix transposition (for details see turtorial part 2).

Since multiple sclerosis is a dynamic disease of the myelin sheaths, the connection weights in the connectome also change dynamically. In the tutorial part 2 we explain how the coupling matrix L was temporally modulated to realize the dynamics of demyelination and remyelination within a simulation: (4) Following the implementing the basic properties of the connectome diffusion models, two methods for generating additive noise and stochastic noise (Ornstein-Uhlenbeck process) were implemented to verify the stability of the reactions-diffusion model with respect to different noise levels. To realize this, a coupled Ornstein-Uhlenbeck process δ was applied additively to the directed diffusion. The different steps for implementing this function were also explained in detail in Tutorial Part 2: (5) The advantage of our connectome frameworks is that neuronal connectivity is embedded in a standard stereotaxic coordinate system, so that a centroid can be estimated from each region and linear Euclidean distances can be estimated between centroids. Of course, real axonal connections run with curvatures as well as partially stronger unsmooth kinks. Therefore we speak here also only of estimations. We have presented in tutorial part 2 an approach to map Euclidean distances to graph edges by inserting nodes. Here, node E has been inserted to specify the distance between node A and node B. D is the diffusion constant and Θ controls the levarge and delay of the inserted node respectively: (6) (7) (8)

We validated the Gray-Scott implementation by applying it to a regular grid with the parameters published by Buric [89] and were able to generate concentration patterns (Fig P in S1 Text) typical of the Gray-Scott model. Furthermore, the Gierer-Meinhardt model was applied to a regular grid and tested with the parameters published in Koch and Meinhardt [56, 58]. As with the Gray-Scott model, the Turing pattern could be generated (Fig Q in S1 Text), which is also known from the literature.

In addition, we tested the runtime of the Gierer-Meinhardt model on a i7–6500U CPU (2,5 GHz). Here, a linear relationship between the number of iterations of the diffusion model and the computation time is found. 106 iterations need about 60 s and 104 iterations about 2 s (the multiple sclerosis mechanosensitive connectome was used). The number of meaningful iterations depends on the parameters of the model, the step size and the structure of the network and must be tested accordingly. Copying the concentration states at the time points for further analysis or graphical representation also takes time, but this was not considered in the values given above. In addition to the influence of the number of iterations, we studied the effect of the number of nodes. Here we found a quadratic relationship between the number of nodes and the computation time for the Gierer-Meinhardt model.

We are pleased to provide the Java source code on request and have also made it available for download at figshare (see below). A direct download is available for the connectivity data and for the executable Java version for Linux, Windows and iOs platforms on the neuroVIISAS web page.

Results

GM and MM in an embedded pathway of a subconnectome

The central pathways of mechanosensitivity originate from the central processes of pseudounipolar neurons in the dorsal root ganglia outside the spinal cord. The connections of the first cervical segments of this pathways were filtered from the complete connectome and both sides of the central nervous system were selected to compute the weighted, directed and bilateral adjacency matrix (Fig 1). These segments manifest particularly severe changes in multiple sclerosis [90, 91]. A symmetric graph representation of this matrix is shown in Fig 2. The source of the mechanosensitive pathway starts in the peripheral nervous system from the mechanoreceptors of the subepidermal layers of the skin or from internal organs. In multiple sclerosis, the mechanosensitive projection exhibits particularly pronounced alterations at specific sites in the spinal cord, primarily causing common neurological symptoms [92]. We consider the central projection from the dorsal root ganglion (first neuron) to the ipsilateral cuneate nucleus (second neuron) then to the right ventrolateral thalamic nucleus (third neuron) with termination in the right somatosensory cortex (terminal neuron). This mechanosensitive pathway is embedded in the bilateral network shown in Fig 2. The first three cervical spinal cord segments were used in the RD modeling.

thumbnail
Download:
Fig 1. The weighted adjacency matrix of a bilateral mechanosensory subconnectome.

The weighted adjacency matrix of the spinal cord, brainstem, diencephalic and cortical connectivity of the mechanosensory pathways. The last character of the area abbreviation indicates the side of the hemisphere: L: left hemisphere, R: right hemisphere. AGl: Lateral agranular prefrontal cortex, AGm: Medial agranular prefrontal cortex, CERC: Cerebellar cortex, Cu: Cuneate nucleus, DCeN: Cerebellar nuclei, DRGC1: Dorsal root ganglion of cervical segment 1, DRGC2: Dorsal root ganglion of cervical segment 2, DRGC3: Dorsal root ganglion of cervical segment 3, Gr: Gracile nucleus principal part, ILN: Intralaminar nuclei, IO: Inferior olive, LTNG: Lateral thalamic nuclear group, mPFC: Medial prefrontal cortex, Pn: Pontine nuclei, PTG: Posterior group, S1: Primary somatosensory cortex, S2: Secondary somatosensory cortex, VL: Ventrolateral thalamic nucleus, VNT: Ventral thalamus, VPL: Ventral posterolateral thalamic nucleus.

https://doi.org/10.1371/journal.pcbi.1010507.g001

thumbnail
Download:
Fig 2. The weighted digraph of the bilateral mechanosensory subconnectome.

The bilateral weighted digraph of the adjacency matrix shown in Fig 1. The dashed lines indicate contralateral projections.

https://doi.org/10.1371/journal.pcbi.1010507.g002

Change of oscillations following reduction of connection weights.

In the following, it was investigated how a single reduction of a connection weight affects the transmission in a coupled RD system. The weight reduction was realized by reducing normalized logarithmic transformed connection weights by 90% and keeping this reduction constant until the end of the process. To perform weighted network analysis, the ordinal weighted categories x were transformed to an exponential scale [93]: (9)

In particular the oscillation of the target region right somatosensory cortex was compared with the unchanged connections weights. The functions of concentrations of the GM RD model are displayed for all regions of the left side (Fig 3) and right side (Fig BJ in S1 Text) in distinct diagrams. Both presentations of concentrations are from the same GM RD-simulation. Right hemispheric (contralateral to initial conditions) regions show lower concentration amplitudes (first amplitudes are < 2.5) than regions on the ipsilateral side of initial condition for the DRGs.

thumbnail
Download:
Fig 3. Functions of concentrations for regions of left side.

The functions of concentrations of left hemispheric regions are shown for the GM RD-model. Initial conditions V0 = 1 and W0 = 1 were set for the dorsal root ganglia of the left side. The color coding of the functions in this diagram and all following function representations in diagram form are based on the color definition in the area hierarchy, which is simplified in the adjacency matrix in Fig 1. Thus, the colors of the columns and rows of the adjacency matrix represent the color scale of the functions. The x-axis shows the iteration steps of the function. On the y-axis the diffused concentrations are shown according to the applied models and functions. These axis assignments were also maintained uniformly for all subsequent function diagrams.

https://doi.org/10.1371/journal.pcbi.1010507.g003

The functions of concentrations in Fig 3 and Fig BJ in S1 Text show variable amplitudes. The largest amplitude (blue) belongs to the cerebellar cortex which does not have direct efferents to non-cerebellar regions. The concentrations of the GM RD are relatively strong due to local connectivity in the cerebellum. Largest logarithmic correlations were found between the coefficient (convergent-divergent degree coefficient) (c = 0.726) (Fig 4 top) and the cluster-coefficient for output connections (CluCout) with c = 0.728. The ranks of local parameters show a relatively low average rank for the cerebral cortex region (Fig 4 bottom).

thumbnail
Download:
Fig 4. Relation of average concentrations and local network parameters.

a) The correlation c = 0.726 of average concentrations (AvgGEMEx) and the CluCOut () coefficient. b) Each point indicates the rank of a local network parameter of a region. The green curve represents the mean ranks. The colors correspond to the color assignments to the considered regions of the network, which were introduced in Fig 1 in the form of an adjacency matrix.

https://doi.org/10.1371/journal.pcbi.1010507.g004

To apply the models to a connectome composed of regions affected in multiple sclerosis disease, the definition of the origins of the somatosensory pathway are particularly important, as this is where the initialization of the models as well as the change in pathway properties is set. The literature shows that the upper cervical spinal cord segments are particularly affected by demyelination [90, 91]. The central axons of the pseudounipolar neurons in the cervical spinal ganglia (1st neuron) project with a branch in the posterior funiculus to the 2nd neuron in the cuneate nucleus and a colateral via a switch in the Rexed layers V-VII (nucleus proprius of spinal cord) via the ventral spinothalamic tract to the ventral posterolateral nucleus of the thalamus. For simplicity, we consider only the projection of the pseudounipolar neurons via the posterior funiculus and, in particular, the cuneate funiculus to the cuneate nucleus. The regions of interest (ROI) are the DRG’s of the first three cervical segments of the left side, which were initialized with V0 = 1 and W0 = 1. The three concentration functions of the left DRGs start at 1 and were documented graphically (Fig BK in S1 Text). In this figure, as previously described, the individual regions were assigned the same colors as in other diagrams with function curves: The small blue curve indicates the concentration in the cuneate nucleus. The magenta curve indicates concentration in the contralateral ventrolateral thalamic nucleus. The brown concentration curve shows the concentration of activators in the contralateral (right hemispheric) somatosensory cortex. The weights of all three connections from left DRG’s 1–3 to the left cuneate nucleus were reduced by a factor of 0.3 for the whole simulation time. It turns out that amplitudes of the ipsilateral cuneate nucleus, contralateral ventrolateral thalamic nucleus and contralateral somatosensory cortex are larger than without weight reduction (Fig BL in S1 Text). Moreover, the oscillation pattern appears to be more regular, especially for the cuneate nucleus. Furthermore, the amplitudes of the contralateral somatosensory cortex increase obviously stronger than those of the cuneate nucleus or ventrolateral thalamic nucleus.

Stable coherent oscillations can be observed even for parameter variations.

If pairs of regions have similar courses of concentrations over a range of different parameters of a RD-model then they might constitute coherent dynamic or functional groups. This has been investigated by varying RD-parameters of the GM and MM models. Here, the cross-correlations of all pairs of regions of the mechanosensitive subconnectome as introduced in Fig 1 has been determined for each variation of parameters. The average cross-correlations over all variations of parameters for each pair of regions is displayed as an average cross-correlation matrix. For the GM model the parameters reaction constants for the reaction substances A (rateA) and I (rateI) as well as the decay rates for the substance A (muA) and I (muI) were investigated over a range from 0.004 to 0.016 with a step size of 0.004. 256 simulations were performed (Fig 5). The average cross-correlation matrix for the variation of the two reaction constants is shown in Fig 6. The average cross-correlation matrix for the variation of the four parameters (2 reaction constants and the 2 decay parameters) is shown in Fig 7. Spectral clustering was performed to detect those regions which share similar average cross-correlations. The group of the 2 parameter reaction constant simulation with 16 simulations consists of left and right hemispheric primary and secondary somatosensory cortex, primary and secondary motor cortex and ventrolateral thalamic nuclei. Same regions build a cluster when the 4 parameters reaction constant and decay rates were varied.

thumbnail
Download:
Fig 5. Average cross-correlation matrix of four parameters of the GM model.

The average cross-correlation matrix of the variation of the four reaction constants is shown. This matrix has been used for further cluster analysis.

https://doi.org/10.1371/journal.pcbi.1010507.g005

thumbnail
Download:
Fig 6. Average cross-correlation matrix of two reaction constants.

The average cross-correlation of the variation of the two reaction constants has been analyzed by spectral clustering. A coherent group of regions (green rectangle) could be determined.

https://doi.org/10.1371/journal.pcbi.1010507.g006

thumbnail
Download:
Fig 7. Average cross-correlation matrix of four parameters.

The average cross-correlation of the variation of the two reaction constants and the two decay parameters for activators and inhibitors has been analyzed by spectral clustering. A similar group of regions (green rectangle) could be determined as found with a two parameter variation.

https://doi.org/10.1371/journal.pcbi.1010507.g007

The 4 reaction parameters A, B, C, D of the MM-model were varied (A: 10 − 15 with step-size: 1, B: 13 − 18 with step size: 1, C: 8 − 10 with step size: 1, D: 0.2 − 0.6 with step size: 0.1) over 540 simulations. The average cross-correlation matrix is shown in Fig 8. Spectral clustering of the average cross-correlation matrix determined the same cluster of coherently activated regions of the left and right hemisphere: S1, S2, AGl, AGm and VL (Fig 9).

thumbnail
Download:
Fig 8. Average cross-correlation matrix of four parameters of the MM model.

The average cross-correlation of the variation of four reaction parameters is shown. Large positive correlation values close to 1 indicate a strong similarity of concentrations of a pair of regions.

https://doi.org/10.1371/journal.pcbi.1010507.g008

thumbnail
Download:
Fig 9. Spectral cluster analysis of the average cross-correlation matrix of the variation of four parameters of the MM model.

The average cross-correlation of the variation of the four reaction parameters has been analyzed by spectral clustering. A similar group of regions (green rectangle) could be determined as found with the GM parameter variation.

https://doi.org/10.1371/journal.pcbi.1010507.g009

The modularity analysis of the adjacency matrix is shown in Fig 10. The regions are distributed over three modules. Cortical left and right hemispheric regions are separated in two different modules. These two modules contain again primary and secondary somatosensory cortex, primary and secondary motor cortex and ventrolateral thalamic nuclei. In addition the VPL, ILN, PTG and LTNG are assigned to these modules. The spectral cluster analysis of the connectivity matching matrix for input and output connections display same regions within a cluster that has been highlighted in Fig 10. Therefore, the clustering of the connectional structure of regions matches the dynamics of coherent regions with synchronous co-concentrations of the reaction-diffusion models. Furthermore, the coherency appears to be stable within the investigated space of parameters.

thumbnail
Download:
Fig 10. Modularity and spectral cluster analysis.

The modularity analysis of connectivity similarity among regions generates 3 modules. The spectral cluster analysis of the connectivity matching matrix reveal a cluster with similar regions like those of the clustered cross-correlation matrices.

https://doi.org/10.1371/journal.pcbi.1010507.g010

Cluster analyses generate groupings of regions with very similar dynamic properties even under the condition of parameter variations.

To prove the effects of changes of parameters on the similarity of the oscillating functions of regions, the variation of parameters are analyzed by spectral [94, 95] and Markov clustering [96]. Spectral clustering techniques are using the spectrum (eigenvalues) of a similarity matrix and the Markov cluster algorithm (MCL) is an unsupervised clustering technique [96]. The following parameters were used (matrix values are normalized between 0 and 1):

  • Maximum value considered zero (≥ 0) [Default: 1.0e − 7]: Test if matrix values are zero. Values below this threshold are set to zero (convergence can be faster)
  • Maximum difference between values (≥ 0) [Default: 1.0e − 7]: If the difference of two values are below the threshold then they are considered as equal (termination of computing)
  • Loop gain (≥ 0) [Default: 0]: Probability for considering a connection on the diagonal. If values > 0 then connectivity loops or self-connections are used. If 0.5 is set as a loop gain, then each value on the diagonal is set to 0.5
  • Inflation exponent (Default: 2.0): Exponent is iteratively applied to matrix elements

In the first step the similarity matrix of either cross-correlations, co-activation or Kuramoto indices are calculated for a simulation. The similarity matrix allows the comparison of cross-correlations, co-activations or the Kuramoto indices between all pairs of regions. Either the cross-correlations, co-activation or Kuramoto index matrices are used to cluster the regions by spectral and Markov methods.

The parameters of the GM or MM models are modified, followed by the calculation of the similarities and the clustering of the regions. To compare the different clusterings the the Jaccard coefficient is used. For every cluster X in clustering A and every cluster Y of clustering B the Jaccard coefficient J(X, Y) = |X ∩ Y|/|X ∪ Y| is calculated. The corresponding cluster for X is the cluster Y that maximizes J(X, Y). If Y is the corresponding cluster for X, X does not have to be the corresponding cluster for Y. As a measurement of similarity of the two clusterings A and B we use the average Jaccard coefficient over all clusters of A and B and their corresponding cluster: (10)

Following parameter variations the similarities of the calculated clusters are visualized in a matrix (Figs 11 and 12).

thumbnail
Download:
Fig 11. Similarities of clusterings referring GM parameter variation.

The reaction parameters and decay rates for substances A and B were varied 256 times. a) RD functions of all regions using default parameters. b) Similarities of the clusterings using spectral clustering (Light: high similarity, dark low similarity). c) Similarities of the clusterings using Markov clustering (Light: high similarity, dark low similarity). d) Magnification of upper left corner of the matrix in b). e) Magnification of upper left corner of the matrix in c).

https://doi.org/10.1371/journal.pcbi.1010507.g011

thumbnail
Download:
Fig 12. Similarities of clusterings referring MM parameter variation.

The reaction parameters A-D were varied 375 times. a) RD functions of all regions using default parameters. b) Similarities of the clusterings using spectral clustering (Light: high similarity, dark low similarity). c) Similarities of the clusterings using Markov clustering (Light: high similarity, dark low similarity). d) Magnification of upper left corner of the matrix in b). e) Magnification of upper left corner of the matrix in c).

https://doi.org/10.1371/journal.pcbi.1010507.g012

In the case of the GM model the variation of reaction parameters rateA, rateI, muA and muI range from 0.004 to 0.016 with a step size of 0.004 leads to 256 combinations of these parameters. For each of the 256 combinations the similarity of functions between pairs of regions has been calculated by the cross-correlation. Referring to this similarity the regions were clustered and the comparison of these 256 clusterings is visualized in a matrix. This matrix shows the pairwise similarity of region clusterings that result from the simulation using a particular selection of reaction parameters.

The RD functions of all regions using default parameters are shown in Fig 11a. The similarities of the clusters using the spectral clustering is shown in Fig 11b. The magnification (Fig 11d) shows a pattern in this matrix. For example the group of parameter sets 5–16 lead to very similar region clusterings. The similarities of the Markov clusterings based on the same parameter variation is shown in Fig 11c. Again the magnification displays homogeneous groups of parameter sets indicating a larger stability of RD functions.

In a comparable way, the similarity, respectively, dissimilarity of clusters of cross-correlations of 375 sets of parameters of the MM model were analyzed. The reaction parameters A, B, and C range from 14 to 18 with a step size of 1. Reaction parameter D ranges from 0.2 to 0.6 with a step size of 0.2 (dP = 0.1 and dQ = 0.01). The RD functions of all regions using default parameters are shown in Fig 12a. The similarities of the clusters using the spectral clustering is shown in Fig 12b. The magnification (Fig 12d) also shows a pattern. For example the group of parameter sets 5–17 lead to very similar clusterings. The similarities of the Markov clusterings based on the same parameter variation is visualized in Fig 12c. The magnification displays highly homogeneous groups of parameter sets indicating a very large stability of RD functions.

Functionally similar regions can be reconstructed based on their dynamic properties in both the GM and MM models.

By comparing the coherent dynamics between MM and GM processes, we want to find out to what extent these two models produce similar results in the same connectome. If similar dynamical behavior of functionally similar regions is obtained in different models, this indicates the reproducibility of a result by another model and thus a certain stability of the dynamics independent of a specific model. Furthermore, we can consider to what extent the similar results of both models can be explained by the specific connectivity rather than by minor changes in parameter settings of the models. We expect that functionally intensive or densely interconnected areas of somatosensory and somatomotor cortical areas will exhibit similar diffusion dynamics and therefore form greater synchronization of concentrations outputs of diffusion functions. This can be tested by comparing pairs of regions with large coherence using a cross-correlation analysis.

The somatosensory regions constitute a set of nodes with common dynamical behavior. A GM and a MM process was applied to the mechanosensory subconnectome. The cross-correlation matrix of co-activations has been determined and analyzed by spectral clustering to obtain groups of regions that share similar synchronous behavior. The left and right hemispheric primary and secondary somatosensory regions build such a group of regions in both RD models (Figs 13 and 14). Moreover, the MM model allows a separation of primary and secondary motoric regions as well as primary and secondary somatosensory regions (Fig 14).

thumbnail
Download:
Fig 13. Result of clustering the cross-correlation matrix of a GM process.

The somatosensory regions constitute a cluster (green rectangle).

https://doi.org/10.1371/journal.pcbi.1010507.g013

thumbnail
Download:
Fig 14. Result of clustering the cross-correlation matrix of a MM process.

The somatosensory regions are contained in the same cluster like in the GM model.

https://doi.org/10.1371/journal.pcbi.1010507.g014

The bifurcation analysis suggests a stable oscillatory behavior of the GM model.

A comparable approach has been suggested elsewhere in order to predict the morphology of pattern generation in a regular lattice [97]. In order to find the oscillation state and the fix-point state of the GM model a simplification using default parameters has been performed. We consider a single PDE without coupling and σ. (11) (12)

Analyzing the fixed points of the GM model by examing the eigenvalues of the Jacobian, reveals the (a, b) plane is divided into an oscillatory region and a region with stable fixed points (Fig 15). For pairs of parameters (a, b) on the right side of the curve r is negative, what indicates stability in the fix-point. E.g., the fix-point u = 2 and v = 4 for the parameters a = b = 0.01 and is stable. This is demonstrated for a single, non connected node in neuroVIISAS (Fig 16a and 16b). The real part r of the first eigenvalue as a function of parameters a and b is displayed in Fig 17. For parameters (a, b) on the left side of the curve the GM model has oscillatory behavior (Fig 18a and 18b). The same parameters were applied to a network with coupling of PDEs. It is demonstrated that the oscillatory behavior is still stable. This becomes apparent from the figures of the functions u and v (Fig BM and Fig BN in S1 Text).

thumbnail
Download:
Fig 15. Linear stability analysis of the GM model.

The curve shows where the real part of the first eigenvalue of the Jacobian matrix in the fix-point is zero in dependence of the parameters a (sigmaa) and b (mua).

https://doi.org/10.1371/journal.pcbi.1010507.g015

thumbnail
Download:
Fig 16. Stability of function u and v.

a) The GM model was applied to exactly one non-connected node to show that the function u has a stable progression above a certain number of iterations. The function of u for a = b = 0.01 shows stability at about 100 iterations. b) For the same parameters as given before, the function v is given.

https://doi.org/10.1371/journal.pcbi.1010507.g016

thumbnail
Download:
Fig 17. Real part of first eigenvalue of the GM model.

The real parts of the first eigenvalue of the Jacobian matrix in dependence of the parameters a (sigmaa) and b (mua).

https://doi.org/10.1371/journal.pcbi.1010507.g017

thumbnail
Download:
Fig 18. Oscillations of functions u and v.

a) In addition to the stability of the function u in Fig 16 the regularity of the oscillations shall be shown here. The regularity extends over all performed iterations. b) The regularity of the oscillations of the function v is shown here. As with the function u, the regularity extends over all iterations performed.

https://doi.org/10.1371/journal.pcbi.1010507.g018

The Wilson Cowan model, like the GM and MM models, leads to stable oscillations with coherent patterns.

The neural mass model of Wilson and Cowan (WC) [98–100] was used (“WC-simulation” in neuroVIISAS [101]) to obtain an visual impression of the dynamics (Fig BO in S1 Text) which results from the same network as used for the GM and MM models (Figs 1 and 2). We used the parameters for a limit cycle oscillation of a single Wilson-Cowan oscillator: aE: 1.2, aI: 2.0, cEE: 5.0, cII: 1.0, cIE: 6.0, cEI: 10.0, θE: 2.0, θI: 3.5, η: 20.0, P: 0.25, E0: 0.1, I0: 0.1, time steps: 1000, time step: 0.1. After about 100 steps the WC model gives rise to stable oscillation within the interconnected regions. The similarity of functions is largest within the peaks of oscillations. The Kuramoto index over all regions displays a regular and stable course. The DRG regions without input connections show a small initialization peak and after about 100 steps they show, as expected, a flat line. Principally, the WC model generates strongly regular oscillations which exhibit some differences when comparing with GM and MM models. The MM model produces much more irregular functions of concentrations with an obvious lower synchronized behavior as can be seen in the smaller Kuramoto indices. The GM model generates more regular oscillations with damping of amplitudes.

The GM and MM models produce different dynamics in degree preserving surrogate networks.

The same number of nodes and connections were used by generating Erdös-Rényi (uniform distributed edges) [102], Watts-Strogatz (small-world) [103], Barabasi-Albert (scale-free) [104], Ozik-Hunt-Ott (small-world) [105], rewiring [106], Klemm-Eguiluz (growing scale-free) [107] and multifractal (cluster coefficient) [108] randomized networks to investigate the effects of structural changes of a network to GM and MM models. The GM model is able to produce regular and stable oscillatory functions within nearly all random networks with the exception of the rewiring network (Fig BH in S1 Text). The cross-correlation matrices display large interregional correlations of function similarities. In the case of the rewiring network the Kuramoto index shows strong changes and functions in early stages of the iterations are damped. Interestingly, the Watts-Strogatz model generates a relatively homogeneous network with regard to edge distribution. However, a chessboard like pattern of large and small cross-correlations can be seen in the cross-correlation matrix. The MM model generates much more irregular oscillations in the different random networks. An obvious feature is the slight overlapping low pass oscillation which is missing in the rewiring network. The rewiring network appears to generate more regularity of the oscillating functions. The limit cycles int phase diagram are lying closer together indicating more similarity of waves (Fig BI in S1 Text).

Dynamic weight changes of the GM model lead to a change in the oscillation of the network nodes of the mechanosensitive pathway.

In the following experiment, the weight modulation has been adapted to a progressive relapsing multiple sclerosis disease progression (Fig 19).

thumbnail
Download:
Fig 19. Weight modulation function.

This weight modulation function (damped cosine function) is applied to the GM RD-system. Initial conditions V0 = 1 and W0 = 1 were set for the dorsal root ganglia of the left side.

https://doi.org/10.1371/journal.pcbi.1010507.g019

By comparing the effect of nonlinear weight modulation of the three connections from the left DRG 1–3 to the left cuneate nucleus with the connectome under normal conditions, the concentrations of the GM RD-system show some remarkable differences. The initial amplitudes of concentrations in the ventrolateral thalamic nucleus and somatosensory cortex are similar (Fig BP in S1 Text) when compared with those in Fig BK in S1 Text. The following amplitudes are slightly more irregular than those in Fig BK in S1 Text with unchanged weights. Furthermore, some small phase shifts are visible when comparing the concentration functions of the contralateral ventrolateral nucleus and somatosensory cortex of nonlinear weight modulated simulation in Fig BP in S1 Text with those in Fig BK in S1 Text. When analyzing the average concentrations by differential connectome analysis [53], strongest differences of the control connectome and weight-reduced connectome are localized at the cuneate nucleus, posterior group of the thalamus, gracile nucleus, lateral thalamic nucleus, intralaminar thalamic nuclei and ventral posterolateral thalamic nucleus. Thalamic nuclei appear to be most affected by a reduction of weights of primary neurons along the mechanosensitive pathway.

To study more specifically the effect of weight reduction within the regions of a pathway embedded in the network, we determined the input regions with diffusion outputs to the regions of the mechanosensitive pathway (Fig BD in S1 Text—Fig BG in S1 Text). It turns out that the ventrolateral thalamus receives slightly changed input following a weight reduction of the afferents of the cuneate nucleus. It is conceivable that such an effect could also occur in multiple sclerosis and lead to transmission disturbances of sensory or mechanosensitive signal patterns with resulting neurological symptoms. These changes are visible after sorting the input nodes by the average concentrations of activators (Fig BG in S1 Text). The ranking of the right cuneate nucleus and the right secondary somatosensory cortex changed (Fig BF in S1 Text in comparison with Fig BG in S1 Text). Stronger changes of ranks of regions with regard to the average concentrations are found for afferents of the right primary somatosensory cortex (Fig BF in S1 Text in comparison with Fig BG in S1 Text).

Long term changes of the oscillation pattern of the regions of the mechanosensitive pathway were observed when iterations were extended up to 10000 (2 [timeStep] × 5000 [timeSteps] = 10000) (Fig BQ in S1 Text). The Kuramoto order parameter r [109, 110] was calculated to estimate the extent of synchronization or desynchronization of the oscillations of concentrations of the regions of the mechanosensitive pathway. The weight modulation function (Fig 20) has been applied to the output connections of the first three cervical DRGs to the cuneate nucleus of the left side. The progression of concentrations of the GM RD-model for the ventrolateral thalamic nucleus and the primary somatosensory cortex of the right side and the cuneate nucleus of the left side is shown in Fig BQ in S1 Text. By applying initial values V0 = 1 and W0 = 1 to the left DRG 1–3 regions the oscillations of the regions of the mechanosensitive pathway develop a slow decrease of amplitudes without change of frequencies or phases. Beside cuneate nucleus, ventrolateral thalamus and primary somatosensory cortex the DRG 2 is plotted in Fig BQ in S1 Text, too. DRG 2 shown an obvious decrease of amplitude size due to its lack of inputs. Following the application of the weight modulation function shown in Fig 20 the amplitudes of concentrations of all these regions of the mechanosensitive pathway decrease until iteration 4700 and then they increase. The Kuramoto order parameter indicates a slight decrease of synchronization of the oscillations of concentrations of the three mechanosensitive regions (Fig BR in S1 Text).

thumbnail
Download:
Fig 20. Weight modulation function.

The weight modulation function (damped cosine function) is applied to the GM RD-system. Initial conditions V0 = 1 and W0 = 1 were set for the dorsal root ganglia of the left side. The number of iterations is 10000. The amplitudes were fixed to 10 and the lower value of the weight reduction was set to 0.1.

https://doi.org/10.1371/journal.pcbi.1010507.g020

Dynamic weight changes of the MM model lead to a change in the oscillation of the network nodes of the mechanosensitive pathway.

The MM RD system works on another time scale than the GM RD. Here, 400 time steps with a step size of 0.1 allows a sufficient survey of the progression of the functions of the concentration of ROIs. To allow comparison of GM and MM the same DRGs of the cervical segments 1–3 of the left side were related to constant non-zero initial conditions with P0 = 6 and Q0 = 12 (Fig BS in S1 Text). Initial amplitudes of the ipsilateral cuneate nucleus, contralateral ventrolateral thalamic nucleus and somatosensory region are slightly delayed or shifted to the right on the x-axis in Fig BT in S1 Text. The amplitudes of concentrations of the cuneate nucleus (blue curve) are better visible in Fig BT in S1 Text because the overlap is not as strong as in Fig BS in S1 Text. The amplitudes of the concentration in the weight-reduced (step function) model of the contralateral somatosensory cortex appear to be smaller than in the control connectome without weight reduction.

A comparable (Fig 19) damped cosine function as a weight modulation function (Fig 21) was applied to the MM RD simulation with the same initial conditions like those 19 in S1 Text. The concentration curves appear to have a lower coherency than in the case of a constant weight reduction (Fig BU in S1 Text). This can be seen particularly well for the function course of the cuneate nucleus, which was assigned the color cyan. Before the reduction of the edge weights, the amplitudes of the cuneate nucleus are hidden by a grouping of amplitudes of most other regions in smaller areas or smaller sections on the x-axis. In other words, most of the amplitudes are co-located and co-occur. If the weights are now reduced, the amplitudes of the cyan coded cuneate nucleus move out of the amplitude group, so that the coherence or synchronization of the amplitudes has slightly decreased.

thumbnail
Download:
Fig 21. The weight modulation function.

This weight modulation function (damped cosine function) is applied to the MM RD-system. It has the same scale of time like the MM RD simulation.

https://doi.org/10.1371/journal.pcbi.1010507.g021

The introduction of distance information creates a damping of the oscillations in the GM model.

The Euclidean distance of the major mechanosensitive pathway from the first dorsal root ganglion and cervical spinal cord segment through the ipsilateral cuneate nucleus, contralateral ventrolateral thalamic nucleus and ventroposterolateral thalamic nucleus to the primary somatosensory cortex (Figs 22 and 23) were used in a GM RD model. When ignoring the distance parameter the oscillation of concentrations is shown in Fig BV in S1 Text. By applying Euclidean distance (d) and log-transformed weights through the following transform F1 × d + F2/w + S (F1 = 0.01, F2 = 0.01, , S = 1) to the GM RD process a damping of oscillations was found (Fig BW in S1 Text).

thumbnail
Download:
Fig 22. Adjacency matrix.

Adjacency matrix of mechanosensitive subnetwork.

https://doi.org/10.1371/journal.pcbi.1010507.g022

thumbnail
Download:
Fig 23. Network representation of the adjacency matrix.

The network of the adjacency matrix of Fig 22 is visualized.

https://doi.org/10.1371/journal.pcbi.1010507.g023

Discussion

Three different RD models were applied to directed and weighted networks and connectomes. It turns out that the GS model does not generate prolonged oscillations which are necessary to investigate the propagation of dynamic signals or the spreading of information through a connectome. The GM and MM RD models are able to generate oscillation patterns in directed and weighted networks. The implemented GM model generates in a regular lattice Turing-like patterns (white eye Turing pattern) [111–113]. The oscillations of the concentrations in the studied networks and circuits are basically similar to the oscillations of RD systems documented in the literature [56, 114]. However, it cannot be said with absolute certainty that there is a perfect match between the oscillation behavior of the procedures we implemented and the same procedures documented in the literature. One reason for this is that different solvers can be used to solve the differential equations and initial conditions can vary. The GM and MM models are stable with regard to additive noise and stochastic noise of the Ornstein-Uhlenbeck process (Fig Y in S1 Text). In addition to these candidate models, there are other promising diffusion models (Barkley, Brusselator, Keller-Segel) that can be applied to networks in the sense of network diffusion [115–120]. However, in the context of this work, we focused on RD models that seemed promising in terms of accounting for weights, distances, and diffusion directions. Therefore, the above-mentioned models were applied to the mechanosensitive pathway embedded in a partial connectome consisting of bilateral regions of the brainstem, diencephalon and cerebral cortex. This subnetwork is of particular interest due to the multiple sclerosis demyelination disease. The mechanosensitive pathway is of special interest because most patients are suffering because of the main symptoms of pain, hypoesthesia and paresthesia [121, 122]. Other complexes of main symptoms (urinary, bowel, musculoskeletal, throat, speech, vision, central symptoms) and their associated pathways are not considered here. A control connectome without changes was compared with a connectome reflecting the demyelination disease progression of the dorsal fascicle of the spinal cord which contains mechanosensitive fibers from dorsal root ganglia to the cuneate nucleus for cervical segments.

Using linear stability analysis [123–131] and bifurcation analysis [131–135] we computed the oscillation state and fix-point state of the GM model. Now it is possible to generate directly oscillatory states of the GM model and study the its dynamic behavior in connectomes. The effect of parameter variations such as rateA, rateI, muA and muI on the coherence of functions between all pairs of regions provides information of the stability of the RD models. We observed similar clustering of coherent dynamics of regions when changing parameters. This stability was found for the GM as well for the MM model. When comparing the Wilson-Cowan neural mass model with the MM and GM models it turned out that the Wilson-Cowan model produces similar oscillations like the GM model. If the WC, GM and MM models are applied to the same network, they produce distinct spectra of amplitudes and frequencies. Furthermore, we detected oscillatory stability if the connectivity changes from the empirical biological connectome to a scale-free or small-world null model. Interestingly, the oscillatory behavior changes strongly when using a rewiring null model.

By varying reaction parameters of the GM and MM models it was found that the same regions give rise to dynamically related groups. The primary and secondary somatosensory and motoric cortical regions as well as the ventrolateral thalamic nucleus of both hemispheres were assigned by spectral clustering to the same cluster. This group constitutes a functional cluster for controlling somatomotor (AGl, AGm) behavior in dependence of sensory input (S1, S2) from the ventrolateral thalamic nucleus.

In the case of a weight reduction by a step function obvious changes of oscillations in the target region of the mechanosensory pathway which is the contralateral primary somatosensory cortex were found. Moreover, it is possible to change the connections weights in analogy to a progressive remitted disease progression of the disease (it is one of the four forms of multiple sclerosis summarized as follows: Clinically isolated syndrome (CIS), Relapsing-remitting MS (RRMS), Primary progressive MS (PPMS), Secondary progressive MS (SPMS)) with distinct changes of oscillations in the contralateral somatosensory cortex [66, 67, 136]. In the meantime, there are models to simulate the clinical course of the aforementioned CIS relatively accurately [71, 137, 138]. In the present study, however, we have restricted ourselves to simply describable functions of weight modulation in order not to consider too many different influencing factors. After applying such a periodic weight reduction (Fig AH in S1 Text), increased amplitudes of the contralateral target region (somatosensory cortex) and slight phase shifts were found. The effect can be carefully interpreted as a dynamic change caused by a structural change in the organization of the network’s connection weights. Moreover, the change of dynamics or oscillation of regions within the mechanosensory pathway may shape the demyelinations effects resulting in paresthesias in the contralateral somatosensory cortex. By investigating local effects of activator concentration of afferents of the mechanosensitive pathway, a change of ranks of regions projecting to the ventrolateral thalamic nucleus and the somatosensory cortex was found. This could indicate a change of input priority of afferent regions or a disruption of temporal signal patterns processed by these thalamic and somatosensory regions. Further, models of weight modulation should be considered as the one suggested by [139]. So far, there exist no computational models that shape the effects of demyelination at the level of networks and connectomes. Most modeling research of multiple sclerosis is performed at the cellular level [71, 137, 138, 140–144], statistically [145] or in MRI studies [146, 147]. We have mentioned above the comparison and application of changes in connection weights with the courses of neurological deficits of CIS of multiple sclerosis. A disadvantage of using RD models is the comparability of time scales to the clinical courses of multiple sclerosis. The clinical courses of multiple sclerosis are over years to decades with oscillatory phases that can last for several months. The RD models are represented by the iterations and iteration step sizes of the functions. The relationship between these fundamentally different scales of neurological progression and function iteration of RD models remains to be worked out.

Further options for extending more parameters to RD models are the consideration of delay differential equations [148–150] to better take note of Euclidean distances in combination with weights of a connectome. A number of works have already investigated the introduction of delays in reaction diffusion systems [133–135, 151], which is also applicable to the future extension of the RD models studied here. Moreover, the consideration of region volume estimates could introduce an important parameter for preferred pathways of the transmission of concentrations of RD through a connectome with many alternative pathways [86, 152].

The implementation of GM and MM RD in the neuroVIISAS framework offers the possibility to either use directional and/or weight information of connectivity. Thus, a RD modeling tool has been made available which enables investigators to apply it to non-directed binary adjacency matrices, complex weighted and directed connectivity data. A further development of these RD approaches studied here is to consider distance-dependent delays in diffusion solved by using delay differential equations which allow the control of a time parameter within a simulation.

Those sets of parameters in the GM and MM models which generate oscillations in empirical weighted and directed networks like the mechanosensitive pathway network give rise to oscillations in other empirical networks with similar global network features. However, it needs to be analyzed how the amount of reciprocal connections, modularity and homogeneity (variability of degrees) may influence RD parameters with regard to the stability of oscillations.

By investigating specific pathways of an empirical network which undergo neuropathological changes like pain pathways, mechanosensory pathways and the visual pathway new techniques are required to judge specific contributions of connections in such pathway-subgraphs. This could support the comprehension how embedded pathways in networks may transmit signals through multiple reciprocal connections as well as alternative routes through connectomes [86]. How signals are transmitted in networks has traditionally been analyzed in gene and protein networks by network propagation and network spreading analysis [153–157]. Signal propagation, signal transmission and signal routing analysis is used to investigate the traveling of signals in generalized and neuronal networks. It has been shown that propagation patterns can be measured in humans and that these patterns are remarkably stable [158]. Therefore, it appears promising to further investigate signal propagation through pathways in partial connectomes shaping certain functions which are disrupted in multiple sclerosis. An important area of work arises from the structured analysis of stimulus response patterns during network diffusion. For this purpose it would be useful to implement FFT [159, 160] and wavelet analyses [161, 162] of the oscillations of nodes of the networks. Several methods exist which demonstrate how analysis of signal transmission may work in neuronal networks [163–166]. Transmission of propagating signals can be analyzed in terms of spreading pattern exhibiting cross networks communications as well as a straight forward transmission with minimal interaction through a certain route of a connectome. At least network flow analysis [167] has been applied successfully for the Caenorhabditis elegans connectome. It should be mentioned that a more basic theory exists of how to model signal propagation in networks using linear response theory [168]. These propagation approaches may be applied to better understand the propagation of signals in non-lesioned control connectomes and connectomes with modulation of connections weights shaping demyelination processes in multiple sclerosis as well as loss of neuron populations in neurodegenerative diseases including Parkinson’s, Alzheimer’s, Huntington’s, Batten disease and amyotrophic lateral sclerosis.

Conclusion

The three reaction diffusion models of Gierer-Meinhardt, Mimura-Murray, and Gray-Scott can be adapted in directed and edge weighted networks or connectomes of the nervous system. Consideration of Euclidean distances between nodes of the treated network in the diffusion processes was also realized. It was shown that the reaction diffusion systems retain their stability against noise. The comparison of the results of the reaction diffusion systems developed here with previously published studies, showed a qualitative agreement. The compound weights can be modulated during a diffusion process. Thus, it is possible to apply demyelination and remyelination processes as they occur in certain forms of multiple sclerosis in a reaction diffusion simulation. It was shown what effects periodic modulations of connection weights have on the oscillation of concentrations in network nodes during diffusion processes. The described implementation was done in neuroVIISAS and is directly executable on different operating systems. In summary, it can be stated that the three reaction diffusion models mimic the effects of demyelination through weight changes within the diffusion process and could be a promising tool to predict changes of connectivity of those regions which do not show obvious function changes.

Supporting information

S1 Text. Supporting information file containing tutorials and supplemental Figs A-BW.

The Supporting Information is divided into several sections and is intended to make the use of the discussed RD models in neuroVIISAS comprehensible. For this purpose, a short introduction to neuroVIISAS is given in the first part of a tutorial. In the second part of the tutorial more detailed background information about the RD models is given. In the third part of the tutorial, complementary experiments to the main findings described in the article are presented. In the last section of the Supporting Information, graphical representations are shown to which references are given in the main body of this paper. They are used for comparison to the described experiments and for the consolidation of individual findings. The exemplar connectivity data (MS.brain) that can be loaded directly into neuroVIIAS and the source code (Reaction-diffusion.zip) are available on figshare (https://doi.org/10.6084/m9.figshare.21081418.v1).

https://doi.org/10.1371/journal.pcbi.1010507.s001

(PDF)

Acknowledgments

We greatly appreciate the assistance by Markus Kipp of the University of Rostock who substantially supports this project and contributed many discussions regarding the disease progression of multiple sclerosis.

References

  1. 1. Deco G, Jirsa VK, McIntosh AR. Resting brains never rest: computational insights into potential cognitive architectures. Trends Neurosci. 2013; 36: 268–274. pmid:23561718
  2. 2. Turing AM. The chemical basis of morphogenesis. Phil Trans R Soc Lond B. 1952; 237: 37–72.
  3. 3. Prigogine I, Lefever R. Symmetry breaking instabilities in dissipative systems. J Chem Phys. 1968; 48: 1695–1700.
  4. 4. Castets V, Dulos E, Boissonade J, De Kepper P. Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern. Phys Rev Lett. 1990; 64: 2953–2956. pmid:10041855
  5. 5. Ouyang Q, Swinney HL. Transition from an uniform state to hexagonal and striped Turing patterns. Nature 1991; 352: 610–612.
  6. 6. Murray JD. Mathematical Biology, I. An Introduction, third edition. Interdiscip Appl Math. 2002; Vol. 17, Springer-Verlag, New York.
    • 7. Barrat A, Barthélemy M, Vespignani A. Dynamical processes on complex networks. Cambridge Univ. Press, 2008.
      • 8. Othmer HG, Scriven LE. Instability and dynamic pattern in cellular networks. J Theor Biol. 1971; 32: 507–537. pmid:5571122
      • 9. Othmer HG, Scriven LE. Nonlinear aspects of dynamic pattern in cellular networks. J Theor Biol. 1974; 43: 83–112. pmid:4813541
      • 10. Horsthemke W, Lam K, Moore PK Network topology and Turing instability in small arrays of diffusively coupled reactors. Phys Lett A. 2004; 328: 444–451.
      • 11. Moore PK, Horsthemke W. Localized patterns in homogeneous networks of diffusively coupled reactors. Physica D. 2005; 206: 121–144.
      • 12. McDougal RA, Conte C, Eggleston L, Newton AJH, Galijasevic H. Efficient simulation of 3D reaction-diffusion in models of neurons and networks. Front Neuroinform. 2022; 16: 847108. pmid:35655652
      • 13. Oliveira RF, Kim M, Blackwell KT. Subcellular location of PKA controls striatal plasticity: stochastic simulations in spiny dendrites. PLoS Comput Biol. 2012 8(2): e1002383. pmid:22346744
      • 14. Stiles JR, and Bartol TM. Monte Carlo methods for simulating realistic synaptic microphysiology using MCell. In: Computational Neuroscience: Realistic modeling for experimentalists, 2001, ed. De Schutter E. CRC Press, Boca Raton, 87–127.
        • 15. Kerr R, Bartol TM, Kaminsky B, Dittrich M, Chang JCJ, Baden S et al. Fast Monte Carlo simulation methods for biological reaction-diffusion systems in solution and on surfaces. SIAM J. Sci. Comput. 2008, 30: 3126–3149. pmid:20151023
        • 16. Wagner M, Bartol TM, Sejnowski TJ, Cauwenberghs G. Markov chain abstractions of electrochemical reaction-diffusion in synaptic transmission for neuromorphic computing. Front Neurosci. 2021; 15: 698635. pmid:34912188
        • 17. Li Y, Kahraman O, Haselwandter CA. Stochastic lattice model of synaptic membrane protein domains. Phys Rev E. 2017; 95: 052406. pmid:28618626
        • 18. Blinov ML, Schaff JC, Vasilescu D, Moraru II, Bloom JE, Loew LM. Compartmental and spatial rule-based modeling with Virtual Cell. Biophys J. 2017; 113: 1365–1372. pmid:28978431
        • 19. Schaff J, Fink CC, Slepchenko B, Carson JH, Loew LM. A general computational framework for modeling cellular structure and function. Biophys J. 1997; 73: 1135–1146. pmid:9284281
        • 20. Novak IL, Gao F, Choi YS, Resasco D, Schaff JC, Slepchenko BM. Diffusion on a curved surface coupled to diffusion in the volume: application to cell biology. J Comput Phys. 2007; 226: 1271–1290. pmid:18836520
        • 21. Blasius TL, Reed N, Slepchenko BM, Verhey KJ. Recycling of kinesin-1 motors by diffusion after transport. PLoS One. 2013; 8: e76081. pmid:24098765
        • 22. Brown S-A, Holmes RM, Loew LM. Spatial organization and diffusion in neuronal signaling. In: Computational Systems Neurobiology 2012, Springer Netherlands, 133–161.
          • 23. Blackwell KT. Paired turbulence and light do not produce a supralinear calcium increase in Hermissenda. J Comput Neurosci. 2004; 17: 81–99. pmid:15218355
          • 24. Ness TV, Halnes G, Næss S, Pettersen KH, Einevoll GT. Computing extracellular electric potentials from neuronal simulations. Adv Exp Med Biol. 2022; 1359: 179–199. pmid:35471540
          • 25. Ellingsrud AJ, Boullé N, Farrell PE, Rognes ME. Accurate numerical simulation of electrodiffusion and water movement in brain tissue. Math Med Biol. 2021; 38: 516–551. pmid:34791309
          • 26. Ellingsrud AJ, Dukefoss DB, Enger R, Halnes G, Pettersen K, Rognes ME. Validating a computational framework for ionic electrodiffusion with cortical spreading depression as a case study. eNeuro. 2022; 9: ENEURO.0408-21.2022. pmid:35365505
          • 27. Sætra MJ, Einevoll GT, Halnes G. An electrodiffusive neuron-extracellular-glia model for exploring the genesis of slow potentials in the brain. PLoS Comput Biol. 2021; 17: e1008143. pmid:34270543
          • 28. Spees WM, Sukstanskii AL, Bretthorst GL, Neil JJ, Ackerman JJH. Rat brain global ischemia-induced diffusion changes revisited: biophysical modeling of the water and NAA MR “Diffusion Signal”. Magn Reson Med. 2022. pmid:35452137
          • 29. Juzekaeva E, Gainutdinov A, Mukhtarov M, Khazipov R. Dynamics of the hypoxia-induced tissue edema in the rat barrel cortex in vitro. Front Cell Neurosci. 2018; 12: 502. pmid:30618644
          • 30. Vorísek I, Syková E. Ischemia-induced changes in the extracellular space diffusion parameters, K+, and pH in the developing rat cortex and corpus callosum. J Cereb Blood Flow Metab. 1997; 17: 191–203. pmid:9040499
          • 31. Ichinomiya T. Frequency synchronization in a random oscillator network. Phys Rev E, 2004; 70, 026116.
          • 32. Boccaletti S, Latora V, Moreno Y, Chavez M, Hwang D-U. Complex networks: Structure and dynamics. Phys Rep. 2006; 424: 175–308.
          • 33. Arenas A, Díaz-Guilera A, Kurths J, Moreno Y, Zhou C. Synchronization in complex networks. Phys Rep. 2008; 469: 93–153.
          • 34. Pastor-Satorras R, Vespignani A. Epidemic spreading in scale-free networks. Phys Rev Lett. 2001; 86: 3200–3203. pmid:11290142
          • 35. Colizza V, Pastor-Satorras R, Vespignani A. Reaction-diffusion processes and metapopulation models in heterogeneous networks. Nature Phys. 2007; 3: 276–282.
          • 36. Colizza V, Vespignani A. Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations. J Theor Biol. 2008; 251: 450–467. pmid:18222487
          • 37. Hütt MT, Armbruster D, Lesne A. Predictable topological sensitivity of Turing patterns on graphs. Phys Rev E. 2022; 105: 014304. pmid:35193278
          • 38. Betzel RF, Griffa A, Avena-Koenigsberger A, Goñi J, Thiran J-P, Hagmann P et al. Multi-scale community organization of the human structural connectome and its relationship with resting-state functional connectivity. Netw Sci. 2013; 1: 353–373.
          • 39. Delvenne J-C, Yaliraki SN, Barahona M. Stability of graph communities across time scales. Proc Natl Acad Sci USA. 2010; 107: 12755–12760. pmid:20615936
          • 40. Goñi J, Avena-Koenigsberger A, Velez de Mendizabal N, van den Heuvel MP, Betzel RF, Sporns O. Exploring the morphospace of communication efficiency in complex networks. PLoS One. 2013; 8: e58070. pmid:23505455
          • 41. Mišić B, Betzel RF, Nematzadeh A, Goñi J, Griffa A, Hagmann P et al. Cooperative and competitive spreading dynamics on the human connectome. Neuron. 2015; 86: 1518–1529. pmid:26087168
          • 42. Abdelnour F, Dayan M, Devinsky O, Thesen T, Raj A. Functional brain connectivity is predictable from anatomic network’s Laplacian eigen-structure. Neuroimage. 2018; 172: 728–739. pmid:29454104
          • 43. Nematzadeh A, Ferrara E, Flammini A, Ahn Y-Y. Optimal network modularity for information diffusion. Phys Rev Lett. 2014; 113: 088701. pmid:25192129
          • 44. Abdelnour F, Voss HU, Raj A. Network diffusion accurately models the relationship between structural and functional brain connectivity networks. Neuroimage. 2014; 90: 335–347. pmid:24384152
          • 45. Poudel GR, Dominguez DJF, Verhelst H, Vander Linden C, Deblaere K, Jones DK et al. Network diffusion modeling predicts neurodegeneration in traumatic brain injury. Ann Clin Transl Neurol. 2020; 7: 270–279. pmid:32105414
          • 46. Raj A, Kuceyeski A, Weiner M. A network diffusion model of disease progression in dementia. Neuron 2012; 73: 1204–1215. pmid:22445347
          • 47. Shadi K, Dyer E, Dovrolis C. Multisensory integration in the mouse cortical connectome using a network diffusion model. Netw Neurosci. 2020; 4: 1030–1054. pmid:33195947
          • 48. Raj A, Powell F. Models of network spread and network degeneration in brain disorders. Biol Psychiatry Cogn Neurosci Neuroimaging. 2018; 3: 788–797. pmid:30170711
          • 49. Park BY, Vos de Wael R, Paquola C, Larivière S, Benkarim O, Royer J et al. Signal diffusion along connectome gradients and inter-hub routing differentially contribute to dynamic human brain function. Neuroimage. 2021; 224: 117429. pmid:33038538
          • 50. Avena-Koenigsberger A, Yan X, Kolchinsky A, van den Heuvel MP, Hagmann P, Sporns O. A spectrum of routing strategies for brain networks. PLoS Comput Biol. 2019; 15: e1006833. pmid:30849087
          • 51. Avena-Koenigsberger A, Misic B, Sporns O. Communication dynamics in complex brain networks. Nat Rev Neurosci. 2017; 19: 17–33. pmid:29238085
          • 52. Schmitt O, Eipert P. neuroVIISAS: approaching multiscale simulation of the rat connectome. Neuroinformatics 2012; 10: 243–267. pmid:22350719
          • 53. Schwanke S, Jenssen J, Eipert P, Schmitt O. Towards differential connectomics with NeuroVIISAS. Neuroinformatics. 2019; 17: 163–179. pmid:30014279
          • 54. Mimura M, Murray JD. On a diffusive prey-predator model which exhibits patchiness. J Theor Biol. 1978; 75: 249–262. pmid:745441
          • 55. Nakao H, Mikhailov AS. Turing patterns in network-organized activator-inhibitor systems. Nat Phys. 2010; 6: 544–550.
          • 56. Gierer A, Meinhardt H. A theory of biological pattern formation. Kybernetik. 1972; 12: 30–39. pmid:4663624
          • 57. Gierer A. Generation of biological patterns and form: Some physical, mathematical, and logical aspects. Progr Biophys Molec Biol. 1981; 37: 1–47. pmid:7244249
          • 58. Koch AJ, Meinhardt H. Biological pattern formation: from basic mechanisms to complex strucutres. Rev Mod Phys. 1994; 66: 1481–1507.
          • 59. Pearson JE. Complex patterns in a simple system. Science. 1993; 261: 189–192. pmid:17829274
          • 60. Lee KJ, McCormick WD, Ouyang Q, Swinney HL. Pattern formation by interacting chemical fronts. Science. 1993; 261: 192–194. pmid:17829275
          • 61. Gray P, Scott SK. Autocatalytic reactions in the isothermal continuous stirred tank reactor: isolas and other forms of multistability. Chem Eng Sci. 1983; 38: 29–43.
          • 62. Gray P, Scott SK. Autocatalytic reactions in the isothermal continuous stirred tank reactor: oscillations and instabilities in the system a + 2b → 3b, b → c. Chem Eng Sci. 1984; 39: 1087–1097.
          • 63. Gray P, Scott SK. Sustained oscillations and other exotic patterns of behaviour in isothermal reactions. J Phys Chem 1985; 89: 22–32.
          • 64. Gray P, Scott SK. Chemical oscillations and instabilities. Non-linear chemical kinetics. Clarendon Press, Oxford, 1994.
            • 65. Hausenblas E, Randrianasolo TA, Thalhammer M. Theoretical study and numerical simulation of pattern formation in the deterministic and stochastic Gray-Scott equations. J Comp Appl Math 2019; 364: 112335.
            • 66. Zhou F, Zhuang Y, Gong H, Wang B, Wang X, Chen Q et al. Altered inter-subregion connectivity of the default mode network in relapsing remitting multiple sclerosis: a functional and structural connectivity study. PLoS One. 2014; 9: e101198. pmid:24999807
            • 67. Sbardella E, Tona F, Petsas N, Upadhyay N, Piattella MC, Filippini N et al. Functional connectivity changes and their relationship with clinical disability and white matter integrity in patients with relapsing-remitting multiple sclerosis. Mult Scler. 2015; 21: 1681–1692. pmid:26041799
            • 68. Ropero Peláez FJ, Taniguchi S. The gate theory of pain revisited: modeling different pain conditions with a parsimonious neurocomputational model. Neural Plast. 2016; 2016: 4131395. pmid:27088014
            • 69. Tommasin S, De Giglio L, Ruggieri S, Petsas N, Giannì C, Pozzilli C et al. Relation between functional connectivity and disability in Multiple sclerosis: a non-linear model. J Neurol. 2018; 265: 2881–2892. pmid:30276520
            • 70. Fleischer V, Radetz A, Ciolac D, Muthuraman M, Gonzalez-Escamilla G, Zipp F et al. Graph theoretical framework of brain networks in Multiple sclerosis: A review of concepts. Neuroscience. 2019; 403: 35–53. pmid:29101079
            • 71. Kotelnikova E, Kiani NA, Abad E, Martinez-Lapiscina EH, Andorra M, Zubizarreta I et al. Dynamics and heterogeneity of brain damage in multiple sclerosis. PLoS Comput Biol 2017; 13: e1005757. pmid:29073203
            • 72. He Y, Dagher A, Chen Z, Charil A, Zijdenbos A, Worsley K et al. Impaired small-world efficiency in structural cortical networks in multiple sclerosis associated with white matter lesion load. Brain. 2009; 132: 3366–3379. pmid:19439423
            • 73. Guye M, Bettus G, Bartolomei F, Cozzone PJ. Graph theoretical analysis of structural and functional connectivity MRI in normal and pathological brain networks. MAGMA. 2010; 23: 409–421. pmid:20349109
            • 74. Muthuraman M, Fleischer V, Kolber P, Luessi F, Zipp F, Groppa S. Structural Brain Network Characteristics Can Differentiate CIS from Early RRMS. Front Neurosci. 2016; 10: 14. pmid:26869873
            • 75. Tewarie P, Steenwijk MD, Tijms BM, Daams M, Balk LJ, Stam CJ et al. Disruption of structural and functional networks in long-standing multiple sclerosis. Hum Brain Mapp. 2014; 35: 5946–5961. pmid:25053254
            • 76. Gamboa OL, Tagliazucchi E, von Wegner F, Jurcoane A, Wahl M, Laufs H et al. Working memory performance of early MS patients correlates inversely with modularity increases in resting state functional connectivity networks. Neuroimage. 2014; 94: 385–395. pmid:24361662
            • 77. Bando Y, Takakusaki K, Ito S, Terayama R, Kashiwayanagi M, Yoshida S. Differential changes in axonal conduction following CNS demyelination in two mouse models. Eur J Neurosci. 2008; 28: 1731–1742. pmid:18973589
            • 78. Wu GF, Alvarez E. The immunopathophysiology of multiple sclerosis. Neurol Clin. 2011; 29: 257–278. pmid:21439440
            • 79. Beecham AH, Patsopoulos NA, Xifara DK, Davis MF, Kemppinen A, Cotsapas C et al. Analysis of immune-related loci identifies 48 new susceptibility variants for multiple sclerosis. Nat Genet. 2013; 45: 1353–1360. pmid:24076602
            • 80. Brändle SM, Obermeier B, Senel M, Bruder J, Mentele R, Khademi M et al. Distinct oligoclonal band antibodies in multiple sclerosis recognize ubiquitous self-proteins. Proc Natl Acad Sci U S A. 2016; 113: 7864–7869. pmid:27325759
            • 81. Brennan KM, Galban-Horcajo F, Rinaldi S, O’Leary CP, Goodyear CS, Kalna G et al. Lipid arrays identify myelin-derived lipids and lipid complexes as prominent targets for oligoclonal band antibodies in multiple sclerosis. J Neuroimmunol. 2011; 238: 87–95. pmid:21872346
            • 82. Sukhinin DI, Engel AK, Manger P, Hilgetag CC. Building the Ferretome. Front Neuroinform. 2014; 10: 16.
            • 83. Shanahan M, Bingman VP, Shimizu T, Wild M, Güntürkün O. Large-scale network organization in the avian forebrain: a connectivity matrix and theoretical analysis. Front Comput Neurosci. 2013; 7: 89. pmid:23847525
            • 84. Stephan KE, Kamper L, Bozkurt A, Burns GA, Young MP, Kötter R. Advanced database methodology for the collation of connectivity data on the macaque brain (CoCoMac). Philos Trans R Soc Lond B Biol Sci. 2001; 356: 1159–1186. pmid:11545697
            • 85. Scannell JW, Blakemore C, Young MP. Analysis of connectivity in the cat cerebral cortex. J Neurosci. 1995; 15: 1463–1483. pmid:7869111
            • 86. Schmitt O, Eipert P, Kettlitz R, Lessmann F, Wree A. The connectome of the basal ganglia. Brain Struct Funct. 2016; 221: 753–814. pmid:25432770
            • 87. Burns G. Neural connectivity of the rat: Theory, methods and applications. PhD thesis, University of Oxford, 1997.
              • 88. Download link of the data used in this investigation: https://neuroviisas.med.uni-rostock.de/data.tar.gz
                • 89. Buric F. Pattern formation and chemical evolution in extended Gray-Scott models. Master of Science thesis in Complex Adaptive Systems, Division of Physical Resource Theory, Department of Energy and Environment, Chalmers University of Technology, Gothenburg, Sweden, 2014.
                  • 90. Bussas M, Husseini M, Harabacz L, Pineker V, Grahl S, Pongratz V et al. Multiple sclerosis lesions and atrophy in the spinal cord: Distribution across vertebral levels and correlation with disability. Neuroimage Clin. 2022; 34:103006. pmid:35468568
                  • 91. Andelova M, Uher T, Krasensky J, Sobisek L, Kusova E, Srpova B et al. Additive effect of spinal cord volume, diffuse and focal cord pathology on disability in multiple sclerosis. Front Neurol. 2019; 10:820. pmid:31447759
                  • 92. Smith KJ, McDonald WI. The pathophysiology of multiple sclerosis: the mechanisms underlying the production of symptoms and the natural history of the disease. Philos Trans R Soc Lond B Biol Sci. 1999; 354(1390): 1649–1673. pmid:10603618
                  • 93. Bota M, Sporns O, Swanson LW. Architecture of the cerebral cortical association connectome underlying cognition. Proc Natl Acad Sci U S A. 2015; 112: E2093–E2101. pmid:25848037
                  • 94. von Luxburg U. A tutorial on spectral clustering. Stat Comp. 2007; 17: 1573–1375.
                  • 95. Ng A, Jordan M, Weiss Y. On spectral clustering: analysis and an algorithm Dietterich T, Becker S, Ghahramani Z (eds.); Advances in Neural Information Processing Systems. 2002; 14: 849–856. MIT Press, Cambridge.
                    • 96. Van Dongen S. Graph clustering by flow simulation. PhD Thesis. 2000; University of Utrecht.
                      • 97. Gonpot P, Collet J, Sookia NUH Gierer-Meinhardt model: bifurcation analysis and pattern formation. Trends Appl Sci Res. 2008; 3: 115–128.
                      • 98. Wilson HR, Cowan JD. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys J. 1972; 12: 1–24. pmid:4332108
                      • 99. Wilson HR, Cowan JD. A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik. 1973; 13: 55–80. pmid:4767470
                      • 100. Byrne Á, O’Dea RD, Forrester M, Ross J, Coombes S. Next-generation neural mass and field modeling. J Neurophysiol. 2020; 123: 726–742. pmid:31774370
                      • 101. Daffertshofer A, van Wijk BC. On the influence of amplitude on the connectivity between phases. Front Neuroinform. 2011; 5:6. pmid:21811452
                      • 102. Erdös P, Rényi A. On random graphs. Publicationes Mathematica. 1959; 6: 290–297.
                      • 103. Watts DJ, Strogatz SH. Collective dynamics of “small-world” networks. Nature. 1998; 393: 440–442. pmid:9623998
                      • 104. Barabási A-L, Albert R. Emergence of scaling in random networks. Science. 1999; 286: 509–512. pmid:10521342
                      • 105. Ozik J, Hunt BR, Ott E. Growing networks with geographical attachment preference: Emergence of small worlds. Phys Rev E. 2004; 69: 026108. pmid:14995521
                      • 106. Maier BF. Generalization of the small-world effect on a model approaching the Erdős-Rényi random graph. Sci Rep. 2019; 9: 9268. pmid:31239466
                      • 107. Klemm K, Eguílez VM. Growing scale-free networks with small-world behavior. Phys Rev E. 2002; 65, 57102. pmid:12059755
                      • 108. Palla G, Lovász L, Vicsek T. Multifractal network generator. Proc Natl Acad Sci U S A. 2010; 107: 7640–7645. pmid:20385847
                      • 109. Schröder M, Timme M, Witthaut D. A universal order parameter for synchrony in networks of limit cycle oscillators. Chaos. 2017; 27: 073119. pmid:28764398
                      • 110. Ódor G, Kelling J. Critical synchronization dynamics of the Kuramoto model on connectome and small world graphs. Sci Rep. 2019; 9: 19621. pmid:31873076
                      • 111. Yan-Ning W, Ping-Jian W, Chun-Ju H, Chang-Song L, Zhen-Gang Z. Turing Patterns in a Reaction-Diffusion System, Comm Theo Phys. 2006; 45: 761–764.
                      • 112. Landge AN, Jordan BM, Diego X, Müller P. Pattern formation mechanisms of self-organizing reaction-diffusion systems. Developmental Biology. 2020; 460: 2–11. pmid:32008805
                      • 113. Hu J, Zhu L. Turing pattern analysis of a reaction-diffusion rumor propagation system with time delay in both network and non-network environments. Chaos, Sol Frac. 2021; 153: 111542.
                      • 114. Marcon L, Diego X, Sharpe J, Müller P. High-throughput mathematical analysis identifies Turing networks for patterning with equally diffusing signals. eLife. 2016; 5: e14022 pmid:27058171
                      • 115. Barkley D Barkley model. http://www.scholarpedia.org/article/Barkley_model, Scholarpedia, 2008.
                        • 116. Barkley D. A model for fast computer simulation of waves in excitable media. Physica D: Nonlinear Phenomena 1991; 49: 61–70.
                        • 117. Barkley D. Euclidean symmetry and the dynamics of rotating spiral waves. Phys Rev Lett. 1994 Jan 3;72(1):164–167. pmid:10055592
                        • 118. Dowle M, Mantel RM, Barkley D. Fast simulations of waves in three-dimensional excitable media. Int J Bif Chaos 1997; 7: 2529–2545
                        • 119. Nicolis G, Prigogine I. Self-organization in non-equilibrium systems: From dissipative structures to order through fluctuations. 1977 Wiley, New York.
                          • 120. Keller EF, Segel LA. Model for chemotaxis. J. Theor. Biol. 1971; 30: 225–234. pmid:4926701
                          • 121. Minagar A. Multiple sclerosis. A mechanistic view. Academic Press, Amsterdam, 2016.
                            • 122. Compston A, Confavreux C, Lassmann H, McDonald I, Miller D, Noseworthy J et al. McAlpine’s Multiple sclerosis. Churchill Livingstone, 2006.
                              • 123. Holmes WR. An efficient, non-linear stability analysis for detecting pattern formation in reaction diffusion systems. Bull Math Biol. 2014; 76: 157–183. pmid:24158538
                              • 124. Getz MC, Nirody JA, Rangamani P. Stability analysis in spatial modeling of cell signaling. WIREs Syst Biol Med. 2018; 10: e1395. pmid:28787545
                              • 125. Lakshmikantham V, Leela S, Martynyuk AA. Stability Analysis of Nonlinear Systems. Systems & Contral: Foundations & Application, Birkhäuser, 2015.
                                • 126. Madzvamuse A, Ndakwo H, Barreira R. Stability analysis of reaction-diffusion models on evolving domains: The effects of cross-diffusion. Discrete & Cont Dyn Sys—A. 2016; 36: 2133.
                                • 127. Sayama H. Introduction to the modeling and analysis of complex systems. Open SUNY Textbooks, 2015.
                                  • 128. Nikravesh SKY. Nonlinear systems stability analysis: Lyapunov-based approach. CRC Press, Boca Raton, London, New York, 2013.
                                    • 129. Nirody J, Rangamani P. An introduction to linear stability analysis for deciphering spatial patterns in signaling networks. arXiv, 1501.06640, 2015.
                                      • 130. Sarfaraz W, Madzvamuse A. Classification of parameter spaces for a reaction-diffusion model on stationary domains. Chaos, Solitons & Fractals. 2017; 103: 33–35.
                                      • 131. Theodoropoulos C, Qian Y-H, Kevrekidis IG. Coarse stability and bifurcation analysis using time-steppers: A reaction-diffusion example. Proc Natl Acad Sci. 2000; 97: 9840–9843. pmid:10963656
                                      • 132. Zhao H, Yuan J, Zhang X. Stability and bifurcation analysis of reaction-diffusion neural networks with delays. Neurocomputing. 2015; 147: 280–290.
                                      • 133. Dong T, Xia L. Stability and Hopf bifurcation of a reaction-diffusion neutral neuron system with time delay. Int J Bifurc Chaos. 2017; 27: 1750214.
                                      • 134. Tian X, Xu R. Hopf bifurcation analysis of a reaction-diffusion neural network with time delay in leakage terms and distributed delays. Neural Proc Lett. 2016; 43: 173–193.
                                      • 135. Lv T, Gan Q, Zhu Q. Stability and bifurcation analysis for a class of generalized reaction-diffusion neural networks with time delay. Discrete Dyn Nat Soc. 2016; 2016: 4321358.
                                      • 136. Ghanbari Z, Gharibzadeh S. Synchrony analysis: application in early diagnosis, staging and prognosis of multiple sclerosis. Front Comput Neurosci. 2014; 8: 73. pmid:25100985
                                      • 137. Kannan V, Kiani NA, Piehl F, Tegner J. A minimal unified model of disease trajectories captures hallmarks of multiple sclerosis. Math Biosci. 2017; 289: 1–8. pmid:28365299
                                      • 138. Elettreby MF, Ahmed E. A simple mathematical model for relapsing-remitting multiple sclerosis (RRMS). Med Hypotheses. 2019; 135: 109478. pmid:31733531
                                      • 139. Pappalardo F, Rajput AM, Motta S Computational modeling of brain pathologies: the case of multiple sclerosis. Brief Bioinform. 2018; 19: 318–324. pmid:28011755
                                      • 140. Coggan JS, Bittner S, Stiefel KM, Meuth SG, Prescott SA. Physiological dynamics in demyelinating diseases: unraveling complex relationships through computer modeling. Int J Mol Sci. 2015; 16: 21215–21236. pmid:26370960
                                      • 141. Pennisi M, Rajput AM, Toldo L, Pappalardo F. Agent based modeling of Treg-Teff cross regulation in relapsing-remitting multiple sclerosis. BMC Bioinformatics. 2013; 14 Suppl 16: S9. pmid:24564794
                                      • 142. Mohan TRK, Sen S, Ramanathan M. A computational model for lesion dynamics in multiple sclerosis of the brain. Int J Mod Phys E. 2008; 17: 930–939.
                                      • 143. Zhang W, Wahl LM, Yu P. Modeling and analysis of recurrent autoimmune disease. SIAM J Appl Math. 2014; 74: 1998–2025.
                                      • 144. Chaubey S, Goodwin SJ. A unified frequency domain model to study the effect of demyelination on axonal conduction. Biomed Eng Comput Biol. 2016; 7: 19–24. pmid:27103847
                                      • 145. Pellegrini F, Copetti M, Sormani MP, Bovis F, de Moor C, Debray TP et al. Predicting disability progression in multiple sclerosis: Insights from advanced statistical modeling. Mult Scler. 2019; 5: 1352458519887343. pmid:31686590
                                      • 146. Pagani E, Rocca MA, De Meo E, Horsfield MA, Colombo B, Rodegher M et al. Structural connectivity in Multiple sclerosis and modeling of disconnection. Mult Scler. 2019: 1352458518820759. pmid:30625050
                                      • 147. Tahedl M, Levine SM, Greenlee MW, Weissert R, Schwarzbach JV. Functional connectivity in multiple sclerosis: recent findings and future directions. Front Neurol. 2018; 9: 828. pmid:30364281
                                      • 148. Otto A, Just W, Radons G. Nonlinear dynamics of delay systems: an overview. Philos Trans A Math Phys Eng Sci. 2019; 377: 20180389. pmid:31329061
                                      • 149. Bellen A, Zennaro M. Numerical methods for delay differential equations, numerical mathematics and scientific computation. The Clarendon Press, Oxford University Press, New York, 2003.
                                        • 150. Kuang Y. Delay differential equations with applications in population dynamics. Academic Press, Boston, 1993.
                                          • 151. Rombouts J, Gelens L, Erneux T. Travelling fronts in time-delayed reaction-diffusion systems. Philos Trans A Math Phys Eng Sci. 2019; 377: 20180127. pmid:31329070
                                          • 152. Girard B, Lienard J, Gutierrez CE, Delord B, Doya K. A biologically constrained spiking neural network model of the primate basal ganglia with overlapping pathways exhibits action selection. Eur J Neurosci. 2020. pmid:32564449
                                          • 153. Biran H, Kupiec M, Sharan R. Comparative analysis of normalization methods for network propagation. Front Genet. 2019; 10: 4. pmid:30723490
                                          • 154. Biran H, Almozlino T, Kupiec M, Sharan R. WebPropagate: A web server for network propagation. J Mol Biol. 2018; 430: 2231–2236. pmid:29524510
                                          • 155. Kalaev M, Smoot M, Ideker T, Sharan R. NetworkBLAST: comparative analysis of protein networks. Bioinformatics. 2008; 24: 594–596. pmid:18174180
                                          • 156. Przulj N, Corneil DG, Jurisica I. Modeling interactome: scale-free or geometric? Bioinformatics. 2004; 20: 3508–3515. pmid:15284103
                                          • 157. Memisevic V, Milenkovic T, Przulj N. An integrative approach to modeling biological networks. J Integr Bioinform. 2010; 7: 120. pmid:20375453
                                          • 158. Mitra A, Raichle ME. How networks communicate: propagation patterns in spontaneous brain activity. Philos Trans R Soc Lond B Biol Sci. 2016; 371: 20150546. pmid:27574315
                                          • 159. Emery BA, Hu X, Khanzada S, Kempermann G, Amin H. Rich experience boosts functional connectome and high-dimensional coding in hippocampal network. bioRxiv 2022.
                                          • 160. Kosciessa JQ, Lindenberger U, Garrett DD. Thalamocortical excitability modulation guides human perception under uncertainty. Nat Commun 2021; 12: 2430. pmid:33893294
                                          • 161. Percival DB, Walden AT. Wavelet Methods for Time Series Analysis. Cambridge Univ. Press, Cambridge, UK, 2000
                                            • 162. Brown E, Kass R, Mitra P. Multiple neural spike train data analysis: state-of-the-art and future challenges. Nat Neurosci 2004; 7: 456–461. pmid:15114358
                                            • 163. Zhao J, Qin Y-M, Che Y-Q. Effects of topologies on signal propagation in feedforward networks. Chaos. 2018; 28: 013117. pmid:29390642
                                            • 164. Hong D, Man S. Signal propagation in small-world biological networks with weak noise. J Theor Biol. 2010; 262: 370–380. pmid:19836404
                                            • 165. Cessac B, Sepulchre JA. Stable resonances and signal propagation in a chaotic network of coupled units. Phys Rev E Stat Nonlin Soft Matter Phys. 2004; 5: 056111.
                                            • 166. St-Onge G, Young JG, Hébert-Dufresne L, Dubé LJ. Efficient sampling of spreading processes on complex networks using a composition and rejection algorithm. Comput Phys Commun. 2019; 240: 30–37. pmid:31708586
                                            • 167. Bacik KA, Schaub MT, Beguerisse-Díaz M, Billeh YN, Barahona M. Flow-Based network analysis of the Caenorhabditis elegans connectome. PLoS Comput Biol. 2016; 12: e1005055. pmid:27494178
                                            • 168. Rülle D. Smooth dynamics and new theoretical ideas in nonequilibrium statistical mechanics. J Stat Phys. 1999; 95: 393–468.