1. Introduction

The mammalian brain is composed of billions of neurons each making thousands of connections to each other via synapses. Neurons receive and transmit information via dendrites and axons. Neurons and synapses form an intricate network that is responsible for the information processing capacity of the brain. The underlying anatomical structure of the nervous system has been studied since the beginning of last century [1].

A connectome of an organism is defined as the mapping of the organism’s neuronal connections. It incorporates information on various types of neurons, the number and position of synapses, synaptic weight and the direction of connectivity. Connectomics is the discipline that involves mapping connectomes of different species. Caenorhabditis elegans is currently the only organism whose complete connectome is available [2]. Due to technological constraints, research in connectomics is limited to small neuronal populations and organisms. However considerable progress is underway in the field of functional connectomics, where, rather than anatomical connectivity, activity-related maps are studied [3]. The mathematical analysis of neuronal networks requires anatomical connectivity information of neurons and their extensions. Graph theory is the branch of mathematics used for the analysis of such networks. Neurons and neurite extensions are represented as nodes and edges. Once the network is quantified, it can be subjected to further analysis. Networks are ubiquitous in nature and are generally classified as regular-networks, where ideally each node in the network has same number of edges, or random-networks, where the edges between nodes are defined by a probabilistic distribution [4]. However, most real-world networks, including neuronal networks, are neither regular nor random, but fall in-between. These networks are called small-world networks with shared properties of both regular and random networks [5] [6].

One way to quantify and compare network functionality is through the analysis of synchronization. Synchronization is defined as simultaneous behaviour and has been observed in the physical and biological systems [7]. There are several examples of synchrony in the nervous system [8]. The first documented instance of synchronization in the human brain was in electroencephalography (EEG) [9] activity which was a result of simultaneous activity of millions of neurons. Neuronal synchronization plays an important role during various behavioural states of the organism. For instance a small network of neurons called central pattern generators, found mostly in the spinal cord, produce rhythmic behaviour even in the absence of stimuli and sensory feedback [10]. Rhythmic behaviour have also been observed in brain regions responsible for working memory and long-term memory [11], episodic memory [12], sleep [13], resting states [14], information processing and learning [15], and pattern recognition [16]. Oscillations also shape the development of neuronal circuitry during neuronal embryogenesis [17]. It plays an important role in the shaping of visual pathways [18], cortical circuits [19], embryonic cells and hippocampal circuits [20]. Disruptions in synchronous behaviour are associated with number of diseases. Epilepsy is a classic example of abnormal synchronization and is attributed to a number of factors [21]. Alzeimer’s disease [22], ADHD [23], bipolar disorder [24] are few diseases that are also associated with abnormal synchronization. While much work has focussed on the functional aspects of neural synchrony, the factors responsible for the origin of synchrony are poorly understood. Understanding the mechanisms supporting neuronal synchrony may be essential for the development of treatment options for disorders associated with disruption in synchrony.

In this study, the correlation between anatomical connectivity and function in developing neuronal networks is investigated. Anatomical connectivity information was traced and multisite recordings were simultaneously obtained. Since an in vivo model setup is not possible due to practical constraints, a neuronal culture in vitro model is employed. It has been reported that dissociated neuronal tissue can self-organize into functional units that are capable of producing spontaneous electrical activity [25]. Neurons explanted from new born rat pups are cultured on transparent dishes embedded with microelectrodes. Fluctuations in extracellular potential at multiple locations are measured. After a viable network is established, neurons and neurites are traced manually from high resolution images at different regions of interest from three to nine Day in Vitro (DIV) of development. Multisite recordings were obtained during the same period. By employing certain assumptions, the network is represented as a simple undirected graph. Data from multisite recordings were analysed to find correlations between structure and function.

2. Materials and Methods

2.1. Experimental Setup and Hardware

A microelectrode array (MEA) is an 8-by-8 grid of microelectrodes embedded within a glass substrate that is conducive to neuronal growth. It is used in studies that require monitoring of local field potential and simultaneous imaging. The diameter of each electrode is 30µm with an inter-electrode distance of 200 µm. The recording area is surrounded by a glass well within which neuronal cells are cultured. Local field potentials from electrodes are transferred to the amplification stage of the hardware via contact pads at the sides of the dish (Figure 1). Cortices were explanted from 1-2 day old Sprague Dawley rats by cervical dislocation and placed in 10ml of basal media immediately after extraction. Composition of basal media is as follows; fetal bovine serum (25 ml), horse serum (25 ml), glucose solution in deionized water (125mg/ml, 1.25ml), glutamine in sterile deionized water (20 mg/ml, 1.25 ml), penicillin-streptomycin solution (1.25 ml), ham’s F12 media (98 ml) and basal media eagle (98.25 ml). The protocol was approved by the Institutional Animal Care and USE committee (IACUC) at Louisiana Tech University where the experiments were performed. After removal of meninges, cells were mechanically dissociated by titration with trypsin and ethylenediaminetetraacetic acid (EDTA). The number of viable cells in the extract is determined by a hemocytometer count and the cells required for the appropriate plating cell density is estimated.

Figure 1. Schmatic layout of a MEA recording dish.It is a 8-by-8 grid of TiN electrodes surrounded by a glass-ring within which the cells were cultured. The contact pads are seen at the edge of the dish. Simultaneous imaging and recording is possible using MEAs

Prior to the day of plating, MEAs were treated with 70% ethanol for 20 minutes and rinsed with sterilized water. They were placed in a sterilized hood overnight with UV light on. Three hours prior to plating, dishes were treated with poly-L-Lysine (PLL) for an hour. Later, PLL solution was removed via pipette and the dishes were washed with PBS saline solution. Dishes were then allowed to completely dry under the BSL II hood.

Once the MEA dishes are ready, the extracted cells were transferred to the MEA dish surface. After trying plating densities ranging from five thousand to four hundred thousand cells per dish, we found that a cell density of two hundred thousand cells per dish was optimal for network survivability, maintenance and imaging purposes. After plating the sterilized dishes, they were sealed with perforated teflon membranes (permeable to air) and were incubated at a controlled environment (37℃, 95 humidity and 5% CO₂ ). After seventy two hours of incubation, dishes are removed from the incubator for recording and imaging. The MEA dish was placed in recording chamber maintained at 37℃. The entire chamber is electrically isolated by enclosure within an electrically-grounded Faraday cage. We were able to record from fifteen electrodes that manually selected from the 8-by-8 grid. The reference electrode was grounded. Through previous trial recordings, we observed that frequency spectrum of basal activity was between one to fifty hertz. A sampling frequency of 1000 Hz was chosen and a recording session typically lasted 40-60 minutes. The same sets of electrodes were used for recordings on subsequent days. Raw signals were pre-amplified, filtering (low pass filter with cut-off frequency of 50 Hz), amplified, digitized and stored for further processing. All analysis was performed using numerical analysis software, MATLAB v7.5.

2.2. Graphical Abstraction

A small section of the culture network (~ 800 µm by 400µm) was chosen for imaging purposes. Multiple bright-field microscopic images (40x magnification) of the section were taken and were combined to form a montage. Presence of healthy neuronal populations, without cellular debris obstruction, in the electrode vicinity was an important criterion for choosing the region. Imaging of the same region was performed from third day in vitro (DIV 3) to ninth day in vitro (DIV 9). Neurons were distinguished from glial cells by their irregular shape and neurite extensions. Since individual synapses cannot be seen at this magnification, it is assumed that two neurons have a link if they have intersecting neurites. Extracted anatomical structural information is abstracted mathematically as a graph with assumptions listed below. A graph, G = (V,E), where ‘V’ is the node set and ‘E’ denotes an edge set. An example of a simple undirected graph is shown in Figure 2.

Figure 2. Illustration of a simple undirected graph with nodes and edges. The node set E, is given by E={a,b,c,d,e,f,g,h} and the edge set is given by V={e1,e2,e3,e4,e5,e6,e7,e8,e9}

For graphical abstraction, the following assumptions were applied:

(a) Edges are assumed to be undirectional and unweighted. The direction of the synapse and the type of synapse (inhibitory or excitatory) were not taken into account.

(b) Multiple neurite connections between neurons and self-innervations were excluded.

(c) A cluster of neurons that are closely spaced ( ~ 1 µm of separation) with each other are assumed to be fully connected i.e., each node in the cluster is assumed to have a link with every other node.

(d) Nodes are assumed to be unweighted. Namely, neurons were not discriminated via their shape, size and type.

An adjacency matrix was used to represent the graph. In the matrix, ‘1’ represents an existence of an edge between two nodes and ‘0’ otherwise. After mathematical abstraction, two network parameters of interest, namely cluster index (CI) and characteristic path length (λ), were computed for all DIV. CI of a node is defined as the ratio between the number of connections the neighboring nodes make with each other and the total possible connections that can exist between them. CI of a network is an indication of how tightly a network is connected locally. Mathematically, it is defined as,

(1)

where, i is the node for which the CI is calculated and e_ij represents an edge between two nodes i and j. The average og cluster indexes of all the nodes in the network is the CI of the network.

The path length between any two nodes is the minimum distance required to reach from that node to any other node. The characteristic path length is average of all the path lengths of the nodes. It is mathematically defined as,

(2)

After manually counting neurites and nodes, standard graph models to compare the culture network against were generated. Based on the neuron and neurite count, a regular ring lattice network, a random graph network, and a small-world network model was generated. CI and λ were evaluated for the networks each day (DIV 3 to DIV 9).

2.3. Measures of Synchrony

Multi-site basal activity recording were obtained from regions of the culture where cell density was most optimal for imaging. Daily recordings were obtained from DIV 3 to DIV 9. After DIV 9, the quality of the culture begins to deteriorate. Correlation coefficient and cross correlation analysis were measures that were used to calculate the degree of synchrony between all possible electrode pairs.

Correlation coefficient (CC):

CC between any two random variables X and Y is defined as the ratio of products to products of their standard deviation.

Where σ_X and σ_Y are the standard deviations of X and Y respectively, μ_X and μ_Y are the mean values.

Cross-correlation:

Measure of similarity between two signals as of function of the time lag between them is the cross-correlation between the two signals. It is mathematically given by,

Where, f and g are discrete signals and ‘m’ is the time lag of the sample interval term.

3. Results

3.1. Network Dynamics

Isolated and infrequent spontaneous electrical activity was observed as early as DIV 3 when neurons began to extend neurites to neighbouring neurons. This activity persisted through all days of recording. Small connected clusters began to emerge by DIV 6 (Figure 3). Graph isomorphism of traced networks during different stages in development is shown (Figure 4). Individual neurons that are isolated from a cluster either migrate towards a nearby cluster and establish connectivity with it or perish. Smaller clusters sometimes join together to form larger clusters. By DIV 9, the clusters are tightly connected to each other, facilitating the propagation of electrical activity that spans throughout the network. This stage of development was also characterized by strong correlated activity across more than two electrodes.

Figure 3. Representation of a small-subsection of the electrode grid. It is overlaid with the traced network. The growth of a neuronal network over the embedded is clearly seen. The green circles represent groups of neurons with tight local connectivity. The red lines are traces of neurite extensions between clusters of neurons

Figure 4. Graphical representation of culture networks from DIV 6 to DIV 9. Clusters increased size by absorbing unconnected neurons in the vicinity. Clusters that were formed earlier during development had the tendency to grow into larger than others. By DIV 8 and DIV 9, only a small portion of the dish are unoccupied by neuron.. Networks were drawn using the graph analysis and visualization software Pajek (Batagelj et al, 2003)

After counting the number of neurons and neurites in the region of interest, for each DIVs, we generated three standard models of graphs; a regular, random and small-world graph with the same number of nodes and edges of our cultured neuronal network. Regular graphs (also called Ring-lattice graphs) are highly ordered with each node in the network having exactly the same number of nodes. Random graphs (also called Erdős-Renyi graphs) are highly irregular and the existence of a link between any two nodes is determined a probability. By counting the number of links and dividing it by the total possible links, we estimated the probability of an equivalent random graph and generated it computationally. A small-world network was generated computationally as described here [4]. Next, the clustering coefficient and characteristic path lengths was tabulated for different DIVs (Figure 5 and Figure 6) and was computed for our culture network, regular, random and small-world networks.

Figure 5. The clustering coefficient of the network is compared with different network models during different days in vitro. Neural cultures networks have a much higher local CC in comparison to other models CC of the clusters were much lower

Figure 6. The characteristic path length of the network is compared with different network models during different days in vitro. Neural culture networks have a much lower value. CC of the clusters was comparable to that of random networks

3.2. Functional Results

Spontaneous spiking activity was changes in the extracellular potential due to the activity of single neurons. The emergence of spikes is shown in Figure 7. The increase is spiking phenomenon at this stage of development coincides with an increase in the overall cluster index of the network i.e., the local connectivity.

Figure 7. Spontaneous activity began to emerge as early as DIV 3. Shown above is the activity from a DIV 6 recording. Clusters of neurons are seen around the electrode vicinity. Neurons are often irregular in shape with neurite extensions and astrocytes are rounder with dark edges. By DIV 6, neurons are in the process of forming tightly packed clusters

Correlation coefficient for each electrode pair was calculated for a time period where maximum activity was observed. During initial stages of development, CC was low (~ 0.3). By DIV 8-9 there was sustained correlated activity and it increased over three-fold, spreading across multiple electrodes. It is suggested that small functionally active clusters are responsible for producing spontaneous oscillation. These activities are spread out across a wider recording region (figure 8) and are indicated by increase in correlation coefficients.

Network activity can be mediated synaptically or nonsynaptically [26]. Synaptic mediation is carried out presynaptic neurotransmitter release into the synaptic cleft and activation of postsynaptic ion channels. Nonsynaptic communication can be mediated by neurotransmitter diffusion. Sometimes released neurotransmitters, like glutamate, diffuse into the extracellular media and stimulate other neurons with which it is not connected synaptically. These two mechanisms differ in timescale. Synaptic events usually take place in order of a few milliseconds and non-synaptic mediation in hundreds of milliseconds. Calculation of the time scale of activity across electrode indicates the type of communication that exists. To check for timescale of activity across the electrodes, cross correlation analysis of four most active electrode pairs on DIV 9 were performed (Figure 9). Electrode pair-wise cross correlation is plotted for each electrode. Maximum correlations corresponded to time lags that were ≤ 5ms.

4. Discussion

Metabolic network of E.Coli [27], nervous system of caenorhabditis elegans nematode worm [2], social networks and computer routing networks [6] are some examples of complex networks in the physical and real world. The interactions between the network’s elements gives rise to rich dynamics that determines the behavior of the system. Networks were typically modeled as either regular or random networks until it were discovered that most real world networks are neither regular nor random, but fall in-between. They are classified as small-world networks. Dissociated neurons in culture organize themselves into networks that produce spontaneous and synchronized oscillations [28]. They retain several functional properties of their intact in vivo counterparts [29]. However, the changes in network characteristics that correlate with the emergence of synchronous activity are unclear. In this study, neuronal network evolution during development was quantified and transitions correlating with synchronized activity were investigated.

Network characteristics:

After plating neurons on MEA dishes, neurites began to emerge within 24 hours. At this stage, neurons appear to be evenly distributed across the MEA dish surface without preferential attachment to any particular region. By DIV 3, neurons extended their neurites to neighboring neurons. During the same period, isolated low-frequency (less than two events per minute) spiking activity was observed. Since spiking activity was absent prior to the development of neurite extensions, it is unlikely for a set of isolated neurons to produce endogenous activity. Neuronal growth factors and signaling molecules play an important role in how neurites connect with other neurons [30]. By DIV 5, the emergence of clusters of neurons was observed (figure 4). Clusters are a small network of tightly packed neurons. At this stage of development, spatially isolated neurons and smaller clusters (~ three to four neurons per cluster) either migrated towards larger clusters or perished (Figure 4). After DIV5 (Figure 5), a rise in CI was observed and by DIV9 it was similar to that of regular-networks (exact values are tabulated in Table 1). Similarly, as the network self-organized, a decrease in the characteristic path length (Figure 6) was noticed, and the numbers were similar to that of a random network (Table 2). By DIV 9, the network also began producing synchronized oscillations over a wider area. We hypothesize that for wide-spread synchronous oscillations; the activity must originate at a cluster near an electrode and spread to other regions within a short span of time. A network that has a high CI and low λ may facilitate this hypothesis. The anatomical connectivity of the culture follows a small-world pattern of connectivity [31]. After the initial phase of activity that was characterized by isolated spiking events, the activity became more sustained but with poor correlation (Figure 8, DIV 7). The average correlation coefficient of the network was 0.3247 at DIV 7, followed by a three-fold rise by DIV8 and DIV9 (figure 8, DIV 8-9). During the same period there was a rise in the CI and lowering of λ. This supports the earlier hypothesis that the evolution to a small-world network characteristic makes the process of sustained oscillations easier.

Table 1. The clustering coefficient values of the different models and the culture network is tabulated below. The first sub-column of G (n,e) is the number of nodes and the second sub-column is the number of edge. The column headings are Ring-lattice (RL), small-world(SM), Culture Network (CN), Erdos-Renyi (ER) and Cluster (Cr) Div G (n,e) RL SW CN ER Cr 4 100 100 0 NA NA NA NA 5 230 429 0 0.1 0.08 0.01 NA 6 151 135 0 0.14 0.32 0.014 NA 7 199 557 0.5 0.4 0.66 0.02 NA 8 159 557 0.5 0.5 0.64 0.04 .05 9 208 672 0.5 0.54 0.69 0.03 .31

Table 2. Characteristic path length values for different models are shown Div G (n,e) RL SW CN ER Cr 4 100 100 25.02 NA NA 5.14 NA 5 230 429 30.07 4.01 NA 8.87 NA 6 151 135 37.51 4.48 1.5 5.82 NA 7 199 557 25.13 3.52 2.07 2.9 NA 8 159 557 20.13 4.13 8.51 2.8 2.1 9 208 672 26.25 4.96 4.02 3.14 4.5

Figure 8. Basal activity from four most active electrodes and the ccorrelation coefficient map for all electrode pairs from DIV 7 to DIV 9. The color bar represents the Correlation coefficient value (dimensionless unit). There was an almost 3-fold rise in average CC (0.32 to 0.96) in the four most active electrodes. The amplitude of the oscillations varied from 20 µv to 40µv and lasted for hundreds of seconds. During DIV 8 and DIV 9, there was an increase in the number of tightly connected clusters (figure 4) suggesting that increase in cluster index contribute to an increase level of synchrony in the culture

Cross correlation of synchronized activity:

Electrode pair-wise cross correlation function was evaluated to verify if distant, but synchronizing regions indeed have short path links (as suggested by the numerical values for λ tabulated in table 2). It was done to ascertain the synaptic distance between synchronizing electrodes. The distance between electrode centers is approximately 100 µm. Since near identical activity was observed across each electrode, as shown by the high correlation coefficient during DIV 9 (Figure 8) and a low maximum time lag (within ± 5ms, Figure 9), it is likely that there are short and direct synaptic connections between the electrodes under consideration. This timescale is consistent with the biophysical mechanisms of synaptically transmitted action potentials and hence it is unlikely to have been mediated by neurotransmitter diffusion (which is much slower and in the order of few hundred milliseconds). Also, if there are multiple synapses between these links, then it is likely that the activity is propagated through the fastest path (i.e., the path containing the least number of synapses). This supports our earlier hypothesis that network with small λ and high CI support and sustain synchronized oscillations.

Figure 9. Cross-correlation function of a small sample of signals across various electrodes is depicted. (Top) a segment of DIV 9 basal synchronized activities from the four most active electrodes are plotted. The highlighted signal is the electrode with respect to which the correlation function is calculated. (Bottom) plot electrode pair-wise cross-correlation (CC) of the four signals as a function of time shift .The maximum activity time lag between synchronizing electrodes is < 5 ms indication a strong synaptic-mediation of activity and existence of short path lengths between electrodes

We also distinguished the connectivity of the neurons and the connectivity of clusters of neurons and noticed that the connectivity between clusters of the network was initially low. However, by DIV9, there was an increase in cluster connectivity increased that corresponded to stronger synchrony (correlation coefficient ~= 0.9) across multiple electrodes. The connectivity between clusters increased without significant changes to the connectivity of nodes within each cluster (Figure 4, DIV 9).

The overall clustering coefficients (of the nodes) were close for DIV 8 and DIV 9 (0.64 and 0.69). However, a transition in the network led to 6-fold increase in the clustering coefficient values of the clusters (0.0542 to 0.3163 from DIV 8 to DIV 9). The increase in the cluster robustness may also contribute to an increase in synchronized activity. The transitions in connectivity between were not caused by an increase in the number of neurites, but by the reorganization of both existing neurites and smaller clusters to form links with other nodes. The ability of the cultures to transition without increasing the number of neurites implies that, with all other factors being the same, a network that more neural pathways can generate synchronized activities more robustly.

5. Conclusions

It appears that dissociated neuronal culture self-organize into networks with small-world properties. Such properties may be important for the origin and persistence of synchronized oscillations. Theoretical studies that involve building computational models of neuronal sub-systems can benefit from such anatomical connectivity information. Simultaneous investigation of network dynamics and function could help understand a variety of neurological diseases such as tumor, Alzheimer’s, epilepsy and autism. Diseases can be modeled as disruption in networks and such models would be useful in developing treatments.