Reaction-diffusion models in weighted and directed connectomes

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.

Tutorial part 1: Starting neuroVIISAS and loading the MS.brain data neuroVIISAS should be started through a MS-Windows batch file (run.bat) or a Linux script file (run.sh) which are located in the installation directory. If a 32 bit JRE is installed and a 64 bit neuroVIISAS.jar is started from a desktop item then the error information of the console will be not available and it appears that nothing happens. It is recommended to use a 64 bit JRE for a 64 bit neuroVIISAS.jar. The run.sh contains a starting command which can be adapted to the size of a connectome project and directly entered on the console: / u s r / j a v a / j r e 1 . 8 . 0 6 6 / b i n / j a v a −s p l a s h : Images / l o g o s p l a s h . png −j a r −Xms2000m −Xmx10000m −Xss24M neuroVIISAS . j a r The command consists of the path to the Java Run Time Environment (JRE) and the flags -Xms, -Xmx and -Xss for starting the Java Virtual Machine (JVM). -Xms defines the amount of memory up to the maximum of -Xmx memory for the neuroVIISAS Java process which is relatively large in the example above: 10000Mb. Also, -Xss is relatively large and assigns the stack size with 24 Mb which is important for huge hierarchies of connectome regions. If neuroVIISAS is started without these parameters and a project file with larger sizes of regions and connections is loaded then it is probable that error messages will occur in the console window and data files are not loaded.
After loading neuroVIISAS the main window appears (Fig A). The MS.brain project file can be loaded by clicking on "File" in the main menu bar (Fig A-E) and selecting "Open project". The references.bib can be loaded in the "Settings" menu and "Change project settings" (Clicking on the button "Choose new bibtex file") allows to specify a reference database, however, it is not mandatory for quantitative rather than qualitative connectome analysis.
Following a click on the "Analysis" menu in the main window menu bar the "Advanced connectome analysis" item must be selected. In the bottom row of the window a list of basic matrix, table and analyses representations is displayed in the form of "Tabs" (Fig C-H). After opening the "Advanced connectome analysis" window the "Adjacency Matrix" tab is active, however, an empty analysis window will be displayed ( Fig B).
The hierarchy or a list of regions of a connectome can be opened in several ways. The straightforward method is to repeatedly press the "+" key until the lowest level of the hierarchy of regions of the connectome has been reached (Fig C-A). Then the connections between these regions can be displayed by pressing the "Enter" key or by clicking on the "Refresh" button in the lower right corner of the window (Fig C-H). A default display of connections (Fig C-B) with a navigation window (Fig C-C) appears. In the default display, the number of connections between pairs of regions that were found in the tract tracing research literature is presented. It can be enlarged by scrolling the mouse wheel. In order to change the matrix display, the "Settings" window must be opened by clicking on the "Settings" icon ( Fig C-E) which is an important feature for most other matrix representations as well. The "Settings" window ( Fig C-F) can be scrolled down until the radio button "Average weight / Most frequent weight" appears (beside the "F" in Fig C). Thereafter, the ordinal coded average weights are represented in a green shaded default color palette (Fig D).

Fig C.
A configured connectome in the analysis window. Hierarchically organized regions of a connectome can be extended by pressing repeatedly the "+" key. After reaching a hierarchy level which contains the nodes or regions with connections the "Enter" key must be pressed to computed the adjacency matrix (default view) A) The extended triangle representation of the hierarchy of regions of the connectome. B) The default representation of connections (edge count matrix). C) The navigator window. D) The data information window. E) The settings button must be clicked to open the "Settings" window F. By clicking in the "Settings" window F on "Average weight / ..." the adjacency matrix representation will change. H) is the tabular bar with the "Refresh" button in the lower right corner. The average weight view in the analysis window. After selecting "Average weight / Most frequent weight" from the "Settings" window the adjacency matrix will be directly updated. The mouse wheel enlarges the adjacency matrix.
In order to apply a RD model to the selected regions and connections as indicated in the "Adjacency matrix" tab "Analysis panels" in the menu bar (Fig C-G) and "Simulation" must be clicked (Fig E). Thereafter, a new empty Tab will be generated and automatically opened (Fig F). By clicking on the "Settings" icon in the RD model window (see also Fig C-E) the parameters of the Gierer-Meinhardt model can be defined. After setting the parameters, the "Refresh" button must be clicked to perform the RD process. After finishing the RD process the table of results of the GM simulation is filled with average activation values and spike count values. In the MiniView the coactivation matrix is shown which can be enlarged by clicking on the radio button "Show coactivation matrix". To obtain the display of frequency space of the GM simulation the button "Function" must be clicked. In the window a diagram is shown ( Fig G) were all concentrations of regions are initially represented by dots. By choosing a linear or cubic interpolation the functions are displayed as curves. A single function or a set of functions can be displayed by selecting the regions in the results table. The RD model window with the "Settings" menu. The "Settings" window has been reduced and scrolled down. Some default parameters of the Gierer-Meinhardt model were modified (rateI, muI, dI). Randomized initial values between 0 and 1 (V 0 , W 0 ) were used. The DoPri45 solver was applied and Laplacian L = A − D out . Time steps were set to 2000 and one time step has the size of 2. The connection weights were used by checkmarking "Weighted" and a normalization of the coactivation matrix was applied ("Relative to activation").  In the next step, the weights of the output connections of the left dorsal root ganglia of the cervical segments 1-3 will be reduced. In the "Settings" window where GM model parameters were adopted the button "Weight modulation" must be pressed. A new interval for the change of weights need to be defined by clicking on "Add interval". In the message window "Input" the end of the x-interval (iteration) has to be set (3600) and in the following message window the end of the y-interval must be defined (0.1) ( Fig  L). Now the cosine function can be select from the List-window. As scale factors x=1.0 and y=0.2 (amplitude damping) were used. Finally, 10 oscillations should be generated. The parameters are applied after pressing the "Accept" button. Before, the "Apply to selected regions" button is pressed, the three regions in the table of results of the GM process have to be selected by clicking on the rows and holding the shift key. Thereafter, by pressing the "Apply to selected regions" button the "Choose an option" window appears were the "Output" button should be pressed. Then, the weight modulation function is applied to the three selected regions, only (Fig M). The "Refresh" button will perform the FM process with weight modulation of the three highlighted regions whereas all other regions will keep there original logarithmically transformed weights (if the "Function" window does not appear after pressing the "Function" button working with a Linux OS/KDE it may be hidden behind other windows: Alt+Tab, then scroll through the list of windows). The results of the GM process with weight reduction of the first three left side DRG regions is displayed in the "Function" window ( Fig N).   Tutorial part 2: Reaction-diffusion systems (RD)

Introductory remarks
The phenomenon of diffusion was discovered two centuries ago by Graham and Fick [1,2]. Fick found that diffusion flux is proportional to the negative of the concentration gradient of particles (Fick law ). Diffusion can be observed in physics with regard to particles (particle diffusion), in chemistry and in biology (macromolecules, ions, cell organelles). Alan Turing [3] developed a generative concept of diffusion which helps to understand how stable patterns (Turing stability) evolve in artificial and natural systems [4][5][6]. Soon after Turing's discovery, the concept of diffusion was abstracted and carried forward to the field of epidemics. Diffusion helps to understand and predict the spread of infections in compartmental models of contagion like the SI (Susceptible-Infectious), SIS (Susceptible-Infectious-Susceptible) and SIR (Susceptible-Infectious-Recovered) models [7][8][9]. Proceeding from the diffusion of infections in epidemic research, diffusion has been applied to generate information flow in networks [10] in order to determine how the topology of a network or the pattern of connections between regions affects diffusion.
The diffusion equation can be used in a more technical context like the reduction of noise in biological and medical imaging to preserve the geometry of anatomical details [11]. However, the diffusion process itself is not used in this application to generate a dynamic process in irregular grids such as planar graphs, networks or digraphs. After all, the measurement of diffusion of water along plasma membranes and myelin sheaths of neurits is fundamental for in vivo and non-invasively diffusion tensor imaging (DTI) to visualize non-directed neuronal connections in living brains [12][13][14][15][16][17].
The microenvironment of brain tissue is structured by membranes and myelin or bundles of myelinated nerve fiber surfaces. These boundaries separate spatial compartments or phases through an interface that permits metastability and/or bistability of different ional and specific molecular concentration patterns. In these compartments, signals of neurons are transmitted by action potentials which are generated and propagated by ion fluxes or ion transport processes, respectively [18,19]. Such electro-diffusion processes modeled (electrodiffusive Kirchhoff-Nernst-Planck framework) [18,20,21] at a microscopic level of neurons constitute a precondition for coding of memory, learning, motion control, vision and auditory information processing at the macroscopic level [22][23][24].
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 [46].
Complex self-organized patterns, such as spreading pulses and fronts, stationary dissipative forms, rotating waves, and turbulences, are produced by reaction-diffusion systems [3,[47][48][49]84]. Several reaction-diffusion models have been discovered to function on networks by interacting organisms occupying network nodes and diffusively transferring information across links [50][51][52][53][54]. Chemical reactors and networks of diffusively connected biological cells or areas can both be represented using reaction-diffusion models [51-54]. Self-organization analysis in complex networks is challenging, therefore it has been limited to non-equilibrium pattern development as synchronization [55][56][57] or epidemic spreading [58][59][60]. Changes in the diffusion constants of activators and inhibitors destabilize the uniform state of a system, resulting in the spontaneous formation of rhythmic patterns (Turing patterns) in chemical processes, biological morphogenesis, and ecosystems, according to Turing [3].  [4,70] and a generative coactivation pattern in excitable neural networks [71]. One hypothesis standing behind dynamic modeling in networks is that connectome topology shapes, organizes and constrains dynamic processes [70]. By adopting a one-dimensional reaction-diffusion (RD) system [4] it could be shown that this model can produce stable patterns in coactivation matrices [70]. These promising results lead to further experiments using stochastic versions of two well-established computational models of brain activity. These are the excitable FitzHugh-Nagumo neuron and SER models which reproduce some patterns of empirical functional connectivity as measured in resting-state fMRI experiments [72,73]. The deterministic SER model is a three-state cellular automaton model which is capable of generating excitable dynamics. The FHN model is similar to diffusion models and generates more distinct dynamic patterns in modular networks than the SER and RD models [69].
RD and SER are not built for incorporating further features of a connectome like connection weights, Euclidean distances and volume estimates of regions or nodes. Furthermore, originally RD and SER were not applied to digraphs with reciprocal edges.
Here, the FHN model has been adapted to generate oscillatory dynamics on a weighted digraph.
Until now, reaction-diffusion models with varied reactions terms in the differential equations have allowed for better parametric creation of individual connectome topologies. In this study, the predator-prey and activator-inhibitor models were used.
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. A detailed tutorial is provided as Supplementary Information.
The weights of neural connections in a connectome are used to calculate the strength of connections between areas in reaction-diffusion models. The number of nerve fibers that have included a tract-tracing substance or an estimate of the number of traced nerve fibers are encoded by the weights of connections in tract-tracing investigations. As a result, estimations of connection weights do not represent the mean thickness of myelin sheaths studied using transmission electron microscopy or other methods. Most tract-tracing investigations use estimates of link weights or densities, which are divided into three categories: weak, moderate, and strong. Furthermore, connections can be defined without any classifications or in weight categories such as "weak to moderate" or "very strong". Schwanke et al. [74] provides an overview of all categories, their interpretation, relationships, and comparisons. The strength of a connection that facilitates transmission of an electrochemical signal across parts of the nervous system may be described in terms of a relationship between the weights of connections and their functional mean. More information may be sent if the connection strength is high. The weights of connection can be linked to neuropathological alterations in nerve fibers, such as those seen in Multiple Sclerosis (MS), when demyelination and axon degradation occur. In this situation, the connection weights will be lowered. The ordinal scaled estimates were always logarithmically converted if connection weights were utilized in this investigation.
The predator-prey model of Lotka-Volterra [75-79, 84] is one frequently used differential equation model to investigate interactions between species for understanding mechanisms of pattern formation. Reaction-diffusion systems modeling predator-prey interactions show a rich spatiotemporal dynamic (oscillatory behavior, spatial pattern developments). Some examples are the predator-prey system of the Segel-Jackson model [96,97], the Leslie-Gower model [80] and the diffusive predator-prey system with delay effect [81, 82].
Mimura-Murray presented an expanded version of Lotka-predator-prey Volterra's model [83][84][85], which has several benefits over the original idea. As a result, the Mimura-Murray model (MM) looks to be promising for network diffusion adaptation.
Gierer and Meinhardt Gierer and Meinhardt (GM) devised an activator-inhibitor reaction-diffusion model [86][87][88]. It comprises of a reaction term with activator and inhibitor parameters that may be used to construct various weighted digraph oscillation patterns. As a result, the GM model looks to be a promising contender for network diffusion applications. Finally, the Gray-Scott model (GS) has been modified to network diffusion. The GS is a standard mathematical model for isothermal autocatalytic reaction with another type of reaction component in the differential equation [89][90][91][92][93][94][95].
The goal of this research was to look at reaction-diffusion models in terms of weighted digraphs and diffusion distances. Synthetic randomized directed networks with maintained edges and nodes were created using data from the connectome of the rat nervous system. In order to develop more realistic coupled dynamic models based on empirical data, the reaction-diffusion models of Gierer-Meinhardt [86-88], , and Gray-Scott [89-95] were adapted to weighted digraphs while taking Euclidean distances between nodes into account.
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.
Another goal of this research was to look at the consequences of changing connection weights in the reaction-diffusion process, since this might be a starting point for modeling neurodegenerative disease development [98][99][100][101][102][103]. Particularly those types of neurodegenerative disorder that affect neural connections rather than gray matter or neuronal perikarya. Multiple sclerosis is a demyelinating condition that progresses in a wide range of temporal and geographic patterns [104][105][106][107][108][109][110][111][112][113]. The framework used here allows for the definition of sets of source nodes or regions and sets of target nodes of interest embedded in a large connectome context in order to systematically investigate signal propagation or information diffusion through highly complex connectional architectures with precisely defined changes in connection weights. As a result, our simulation environment looks to be a good starting point for investigating changes in network properties in morphopatholgical neurological disorders and their influence on dynamic pattern alterations in a consistent and repeatable manner. General form of a RD model and the Laplace operator RD systems can be used to model the heterogeneous progression of concentrations of interacting substances on a two-dimensional domain. Instead of concentration the terms substance, amount of substance or abundance can be used. The general form of an n-dimensional reaction-diffusion system is: with Y ∈ R n . The Laplace-operator is which allows modeling of spatial diffusion of the two substances. Equation (14) is a partial differential equation which is transferred into a system of ordinary differential equation by spatial discretization. Here, are reaction terms considered which consist of two substances, respectively, the concentrations u and v. First, the case of a two-dimensional system is considered which will be extended to a reaction-diffusion system of networks. For the two spatial dimensions, the second derivatives are approximated by finite differences It is assumed that the distance between two points on a regular lattice is equidistant, then ∆x = ∆y = h and we obtain: To the nodes at the border, periodical boundary conditions are applied. In an example of a regular lattice of 4 × 4 nodes described by u i (t) = u(x i , y i , t) the following unknown functions are searched for: Lattice with numbering of nodes. This is an example of a possibility to number consecutively the nodes of lattice. It is used in the following.
October 11, 2022 23/63 As an example, for the first upper left node (sub-index: 1) in the lattice (Fig O) we obtain A homogeneous linear differential equation with constant coefficients of first order is used here, for example The differential equation systems can be written by using a matrix A = 10 −5 6 −1 or simply as an adjacency matrix A dy dt = Ay.
By applying the discretization through the Laplace operator, the system of partial differential equations is formed to a system of ordinary differential equations: . . .
The constraint should be kept in mind that the application of spatial discretization allows an approximation of a solution for the nodes and not an exact solution.

Numerics of the differential equations
In the following is a system of ordinary differential equations. By setting x(t 0 ) = x 0 the equation (18) presents an initial value problem. An autonomous differential equation system is obtained if the right side f (x, t) = f (x) becomes independent of time. Numerical integration of ordinary differential equations can be performed by temporal discretization and an approximation of the temporal derivative by finite differences. In the following ∆t is the distance between two nodes on a regular lattice of time points. By calculating the feed-forward difference the temporal derivative can be approximated by Consequently, the explicit Euler-solver can be applied.
More generally an explicit single-step technique is presented by with the function Φ. If Φ contains dependencies of x(t + ∆t), then an implicit method is obtained.

Runge-Kutta method
In the Runge-Kutta method the function Φ has the form Φ = b j · K j whereas K j presents the intermediate steps. For a s-step Runge-Kutta-technique and j = 1..s we get: The coefficients c j , b j and a jl were taken from a Butcher-tableau (RK41) [1].

Embedded Runge-Kutta method
Within each step of time the step width ∆t should be small enough in order to stay under a tolerance tol of precision. The error can be estimated by a procedure of orders p and q with q < p. x 1 is a step within the procedure of p-th order and x 1 a step within the procedure of q-th order. Then the error can be estimated by err = x 1 − x 1 ∆t A step is accepted if the tolerance is larger than the error. In the other case the step width will be reduced. The embedded Runge-Kutta technique is step-controlled like the Dormand-Prince procedure which has been applied.

Gray-Scott model
The Gray-Scott model allows the description of a chemical reaction of the form U + 2V → 3V V → P by using ordinary differential equations. The substance V can be considered as an activator and U is a substrate and does not operate as in inhibitor. By applying the mass action law, the following equations are obtained, where F is a constant term of production and −F u or −kv − F v are linear decay terms: As an example a quadratic domain [0, 1] × [0, 1] is considered. The solution of the PDEs has been realized as described in (14). The following parameters were used: Tab 1. Parameter for the Gray-Scott model. The default parameters are the same as proposed by Buric [2].

Gierer-Meinhardt model
The Gierer-Meinhardt model describes the diffusive interaction of two substances. In the following the two-dimensional case is considered. The concentrations of the two substances are given by U (x, y, t) and V (x, y, t) at point or node (x, y) within time point t. The change of the concentrations within time is described by the two PDEs provides the diffusion or spatial spread based on Brownian motion of the two substances. The model depends on 9 parameters, whereas the parameters r u , r v , µ u , µ v , κ, σ u and σ v determine the local reaction of the two substances and the parameters D u and D v are the velocity of spatial spreading of the two substances. The reaction constants of the activator r u and inhibitor r v have default values of 0.01 and 0.02. The decay rates of activator µ u and inhibitor µ v have default values of 0.02 and 0.01 the decay rates for the activator (u) and inhibitor (v). For controlling the level of saturation or limitation of the autocatalysis of the activator κ is used with a default value of 0.025. The constant production rate of the activator is σ u and the constant production rate of the inhibitor is σ v . Firstly, in analogy to the Gray-Scott model the outcome of RD is considered. By using the parameters Tab 2. Parameters of the local reaction. The parameters are based on Koch and Meinhardt [3,4].
and the diffusion constants a Turing-pattern can be generated (Fig Q).

Mimura-Murray model
The prey-predator model of Lotka and Volterra [5] uses two population sizes of preys P (t) and predators Q(t) where P is the predator and Q the prey population. A generalized model has been introduced by [6,7]. The parameters a, b, c characterize the prey population and parameter d belong to the predator population. Here, the two populations are considered on a domain with two spatial dimensions in such a way that the PDEs have the form [8]: The following parameters turned out to produce an equilibrium state in a regular lattice Tab 4. Parameters of the local reaction of the Mimura-Murray model were used as suggested by Asllani et al. [9].

Reaction-diffusion models in networks
Networks are discrete structures and the Laplace operator was written in the form of a matrix after approximation of the second derivative by finite differences. This matrix corresponds on the non-diagonal elements of the negative adjacency matrix of the regular lattice. In the following, the Laplace operator is modified for the application of diffusion in networks. The Laplace matrix of a non-directed graph G = (V, E) is defined as whereas d i is the degree of node i. For a directed graph [10,11] the following possibilities should be differentiated: These variants have been implemented in neuroVIISAS. The difference between L 1 and L 2 is demonstrated in an empirical network which was used for the following experiments as well (Fig R). To allow a comparison with L 1 and L 3 of output diffusion, the equilibrium state of a Turing dot pattern is shown in Fig AT for   When performing a diffusion process without reaction then a single Euler step with standardized step width and neglection of the diffusion constant is:

Weight modulation within diffusion processes
Neurodegenerative disease in particular Multiple Sclerosis leads to time-dependent progressive patterns of functional loss. This functional loss is shaped by changing connectivity weights within the run time of a reaction-diffusion process. That means that the coupling matrix L is time dependent and the reaction-diffusion with time-dependent coupling is: The weight modulation (Fig S) has been realized for subsets of nodes of a network in order to specify diseased nodes and their input or output connections.
To illustrate the change of connectivity weights a small cyclic 4 node graph (Fig T) should be considered. The node 2 is selected for weight reduction of outgoing or efferent connections.
The corresponding adjacency matrix has the form By analogy reduction of connection weights can be defined for ingoing connections as well as ingoing and outgoing connections.

Additive noise
In order to assess the stability of the reactions-diffusion systems a noise term has been added.

Adding internode Euclidean distances to the diffusion
We have described how weighted connections are incorporated into the diffusion process. Furthermore, the distance between regions as the Euclidean distance between gravitation centers of the regions which are characterized by contours in the digitized stereotaxic atlas is regarded in the diffusion term. Large connections weights transfer concentrations by diffusion strongly. Long distances between nodes delay the transmission by reducing the diffusion velocity of matter. Hence, the coupling should be influenced by distance, only. This can be realized by applying delay differential equations (DDE) [12][13][14][15][16][17] to the preceding concentration. Therefore, it is necessary to initiate the systems with more initial conditions because for each node an initial interval with the length of the delay must be provided.
An alternative has been chosen here. Each edge is subdivided by inserted nodes. In the most simple case, a graph with two nodes and one connection is divided by one inserted node (Fig U). The differential equation for modeling the example graph (Fig U) is given by using the diffusion constant D: The system of differential equations of the graph with one inserted node is: The parameter Θ controls the leverage or the extension of delay of the inserted node. The effect is shown by applying the initial values (A 0 , B 0 ) = (1, 0) and (A 0 , E 0 , B 0 ) = (1, 0, 0) for the inserted node, respectively. The diffusion constants are D = 1 and e 1 = 1 in both cases. In Fig V small Θ lead to a slow increase of concentration in node B because the node E interacts as an inserted node. For a large Θ, the inserted cache is depleted faster. Therefore, with increasing Θ the effect of inserted nodes decreases. In the limit case Θ → ∞ the effect of inserted node generated the same progress of concentration for A and B like in the case without inserted node. Especially for digraphs, reciprocal edges should be treated independently. If the edge between 2 nodes is reciprocal then an inserted node should be generated for each direction. For a larger example (Figs. Fig AB, Fig AC) the procedure is described in the supplement.
Stochastic noise A stochastic noise model was used to investigate the stability of the GM reactions-diffusion system. The coupled Ornstein-Uhlenbeck process is compared with additive normal distributed noise.
Coupled Ornstein-Uhlenbeck processes in networks In the following, coupled Ornstein-Uhlenbeck processes in networks are considered in combination with a diffusive coupling. For a diffusion without reaction term the step U t → U t+∆ t is given by whereas L is the Laplace-or the general coupling matrix of the network. The Euler-Maruyama procedure for the start vector U 0 can be written as:

OU-processes with diffusive coupling
The parameters σ = Θ = 1, D u = 0.1 and µ = 10 were defined. As an initial vector random values between 0 and 1 were used. The resulting progression of concentrations is shown in Fig W. For a better overview, the progression of concentration is displayed for a subset of nodes, only. Gierer-Meinhardt with OU-process as a term of noise First of all, the reaction term of the Gierer-Meinhardt reaction-diffusion systems with adding the noise of an OU-process is considered. To this end the stochastic differential equation system can be considered: In the case of σ = 0, the Euler-Maruyama procedure can be simplified to an Euler-procedure with constant increment. The difference of the components U and R is interesting because R presents the following noise version U . The equation of R is built so that R pursues asymptotically against U (t). Subsequently, a delay between U and R appears, which is shown in the following example.
The example function has the form of In Fig X the progression of R for different Θ is shown. In all cases the tendencies to the value of U are recognizable, however, for small values of Θ large differences are obvious. Intuitively accessible Θ represents a measure for the change of R. If the change of dU dt is too large, then it is not possible anymore that with a small Θ an asymptotic convergence of R against U can be reached. In this case the progression of the function of R depends of U . Thus, the parameter Θ can be interpreted as a delay. If Θ becomes sufficiently large, then the derivation of R to U becomes marginal.
Gierer-Meinhardt and Ornstein-Uhlenbeck with σ = 0 In the following Θ should be large enough so that no significant delay between R and U can be expected. The following system is considered This relates to the matrix notation as follows: The differential equation system of the graph with inserted nodes (Fig AC) has the form of: This system can be simplified with regard to matrix notation by using the incidence matrix for the source graph (Fig AB). The edges are indicated by The adjacency matrix of the new graph with inserted nodes (Fig AC) has the form of: With the related Laplace-matrix L 3 the differential equation system describing the diffusion can be written as follows:

Formalization of a RD with regard of distances
In the formalization of the RD with delay or attenuation of the diffusion term it should be emphasized that the reaction term is applied to the source node, only. Proceeding from a RD system with the adjacency matrix A, reaction term f , Laplace matrix L 3 and vector of concentration Y as a differential equation system we get and are able to determine the adjacency matrix A neu and Laplace-matrix L neu 3 of the new network with inserted nodes. As far as the source network contains n nodes and m edges the differential equation system of the extended network has the form: For the implementation of larger networks, the Laplace matrix for the extended network has the dimension (n + m) × (n + m). By taking into account that newly inserted nodes possess only one input and one output, the Laplace matrix is sparse and computing with vectors is faster.

Tutorial part 3: RD systems in node chains and circuits
The RD systems of Gray-Scott (GS), Gierer-Meinhardt (GM) and Mimura-Murray (MM) show different progressions of their functions in directed networks with connections weights as well as without connection weights. The results of the RD systems under different initial conditions and constraints with regard to non-weighted, weighted, weight changes and inter-node distances in node-chains, node-circuits and a subconnectome will be presented in the following.

Unique initial conditions
All RD systems were applied to a 6 node chain ( Fig AE). The color code of the nodes were applied to the graphs in the diagrams (Figs. Fig AF-Fig AH). The termination of RD processes were done until multiple detectable peaks of oscillating functions were recognizable. For example, we found this pattern in the GS process after 10000, in the GM process after 5000 and in the MM process after 5000 steps. However, the step sizes in the Runge-Kutta45 solver are adapted dynamically. The GS functions of the 6 nodes show only a few waves which reach a stable equilibrium very fast (Fig AF). The GM function follows a damped oscillation (Fig AG) whereas the phases stay relatively constant. The MM functions display sharp oscillating curves without damping (Fig AH). Peaks of the MM functions curves seem to follow a low-frequency oscillatory process.

Initial conditions for selected nodes
When applying initial conditions (min S = 1, max S = 1, min A = 1, max A = 1) to node 0 in the GS model all functions of the nodes have the same progression. A typical first wide wave with the largest peak is generated followed by a strong damping and 2 succeeding peaks (Fig AI). By assigning the initial condition to node 0 in the GM model (V 0 = 1, W 0 = 1) the oscillations as observed in (Fig AG) are damped. The functions show a relatively constant progression with similar phases. The function of node 6 has the largest initial peak and the strongest damping with small amplitudes (Fig AJ). The GM model appears to be more sensitive to an initial condition for a particular node than the GS model. After setting the initial values of the MM model to min X = 1, max X = 1, min Y = 1 and max Y = 1 the oscillation of all nodes (Fig AK) are comparable like in Fig Fig AH. The nodes 1 to 5 which have initial values min X = 0, max X = 0, min Y = 0 and max Y = 0 generate larger initial peaks than node 0 ( Fig  AG). Because the GS model does not process prolonged oscillations in a network it was not investigated under the following conditions and networks (Figs. Fig AF, Fig AI).

Non-weighted circuit
The GM and the MM models were investigated in a 3 node circuit with the feeder node 0 ( Fig AL). Same sets of parameters for RD equation were used as in the chain graph. Here, the feeder 0 was initialized with parameters V 0 = 1 and W 0 = 1, whereas node 1 gets input by diffusion from the feeder node. The short delay of peaks between node 0 and node 1 can be seen in (Figs. Fig AM, Fig AN). The subsequent nodes 2 and 3 show a small delay to node 1.   The MM model shows a more regular oscillation without decrease of amplitudes ( Fig  AP). Node 1 is the only node in the graph which receives two inputs. This node shows a  (Fig AN).  (Fig AQ). The different oscillations are less in phase than those of the GM RD model (Fig AO).
GM RD produces patterns in a 32 × 32 lattice A GM RD process was performed on a 32 × 32 lattice with reciprocally connected 4 neighbors (Fig AD). In a topologically ordered matrix the concentrations of each node show a discrete pattern of high and low activator concentrations (Fig AT). The inhibitor concentrations generate a stable Turing pattern as well (Fig AU). The implemented GM RD functions produce regular dot type Turing patterns that are similar to those shown above when using the same parameters (see parameter lists in Figs. Fig AT and Fig  AU) and sufficient iterations (2000 iterations with 4 iterations: 8000 points of time). The four curves of concentrations are from the 4 corner nodes of the lattice and the center node (orange) (Fig AD).

Weighted circuits and decrease of weights
The weights of connections were applied to the circuit example ( Fig AL) in the GM and MM models. Node 0 has the same initial conditions in the GM and MM models in the following examples. The weighted GM functions (Fig AW) show the damped oscillatory progression which is preserved by introducing connections weights to the RD process. The weighted MM functions (Fig AY) preserve their high-frequency oscillatory progression without damping and relatively constant amplitudes. Then the weight of connection from node 1 to node 2 was reduced by a factor of 10 for the GM and MM models. A relatively strong decrease of amplitudes of node 2 (green) and 3 (blue) can be seen (Fig AX). However, node 1 show brighter waves. A singular reduction of a connection weights obviously affects the pattern of oscillations of this circuit. The MM functions graph of node 1 shows a decrease of amplitude without increase or decrease of frequency (Fig AZ). Functions of node 2 and 3 appear unchanged with regard to the non-reduced connection weight circuit (Fig AZ).

Stability of MM RD with regard to stochastic noise in a subconnectome
The MM RD is relatively stable when applying additive stochastic noise to an empirical network (directed, weighted). The concentration functions are smooth without additive noise (Fig BA). When adding noise with σ = 0.1 the noise effect is obvious at the bottom of the concentration curves (Fig BB). The MM RD functions stay relatively stable when adding noise with a σ = 0.5 (Fig BC) (Fig AL). The weight of the connection from node 1 to node 2 has been reduced.   (Fig AL). The weight of the connection from node 1 to node 2 has been reduced.