Perineuronal nets restrict transport near the neuron surface: A coarse-grained molecular dynamics study

Perineuronal nets (PNNs) are mesh-like extracellular matrix structures that wrap around certain neurons in the central nervous system. They are hypothesized to stabilize memories in the brain and act as a barrier between cell and extracellular space. As a means to study the impact of PNNs on diffusion, the nets were approximated by negatively charged polymer brushes and simulated by coarse-grained molecular dynamics. Diffusion constants of single neutral and single charged particles were obtained in directions parallel and perpendicular to the brush substrate. The results for the neutral particle were compared to different theories of diffusion in a heuristic manner. Diffusion was found to be considerably reduced for brush spacings smaller than 10 nm, with a pronounced anisotropy for dense brushes. The exact dynamics of the chains was found to have a negligible impact on particle diffusion. The resistance of the brush proved small compared to typical values of the membrane resistance of a neuron, indicating that PNNs likely contribute little to the total resistance of an enwrapped neuron.

1. Introduction

Perineuronal nets (PNNs) are mesh-like extracellular matrix structures that ensheath the soma and proximal dendrites of neurons in the brain, particularly parvalbumin-positive (PV) interneurons (van 't Spijker and Kwok, 2017). They have garnered particular interest over the last two decades because of their proposed role in stabilizing long-term memory (Pizzorusso et al., 2002; Tsien, 2013), a hypothesis which have later been corroborated by experiments (Thompson et al., 2018; Carulli et al., 2020). The nets have a complex structure and consist to a large degree of glycosaminoglycans, a family of long, unbranched and highly negatively charged polysaccharides. Due to their large negative charge, the nets are hypothesized to restrict ion transport (van 't Spijker and Kwok, 2017). Developing an improved understanding of how PNNs can act as transport barriers is important because ion concentration and transport affect cell function.

PNNs may be considered to comprise a porous system, with similarities to polymer brushes. However, the exact structure of PNNs remains difficult to characterize. The impact of the structure of PNNs on transport into, out of and around neurons is therefore sparsely understood. Given the limited insight into the precise structure of PNNs, it is important to develop an improved understanding of the general properties of diffusional transport in polymer networks and brushes that reflect general aspects of PNNs. While the fundamental process of diffusion has been well understood since Einstein's seminal paper in 1905 (Einstein, 1905; Tartakovsky and Dentz, 2019), our understanding of diffusion near surfaces and inside tight porous media is still being developed. For example, recent studies address surface diffusion (Postl et al., 2022), diffusion inside nanotubes (Zelenovskiy et al., 2020) and inside planar polymer brushes (Wetzler et al., 2020; Luo et al., 2021) or diffusion of spherical brushes in a polymer melt (Chen et al., 2022). Diffusion in polymer brushes are also of interest in other applications, such as for our understanding of lubricants (Spirin et al., 2011) and for nanoscience and technology (Zhang and Xiang, 2013).

Computational studies of diffusional transport can complement experimental studies in particular in cases where detailed measurements prove difficult, like inside shearing brushes (Ou et al., 2012). Molecular dynamics (MD) simulations have the advantage of yielding deterministic, time-dependent trajectories. A number of molecular dynamics studies have been conducted both on polymers in general and polymer brushes in particular. Singh et al. (2015) applied implicit solvent molecular dynamics simulations through the Langevin equation on apposing neutral polymer brushes represented by a coarse-grained bead-spring model without bending terms in the potential. The density profiles, compression curves and friction was found. Studies of charged apposing brushes using implicit solvent were performed by Ou et al. (2012) and Cao et al. (2010). Both papers found the friction coefficient μ, monomer density profile ρ_m and degree of interpenetration I under shear motion for two apposing charged polymer brushes with bond stiffness and steric effects through the finitely extensible linear elastic (FENE) (Kremer and Grest, 1990) and Lennard-Jones (LJ) potentials, but no potential accounting for bending stiffness. The osmotic pressure and virial terms were plotted against the separation distance by Cao et al. (2010). Ou et al. (2012) found the normal pressure for three different separation distances, but put effort into comparing μ, ρ_m and I of apposing charged and apposing neutral brushes with different counterion valency. Singh et al. (2016) studied the effect of chain stiffness, grafting density and wall separation on apposing electrically neutral polymer brushes by incorporating a tunable cosine bending potential, finding the normal stress, shear stress and friction coefficient. Additionally, they probed the structure by measuring the brush height and number density as a function of chain stiffness and grafting density. A relation between the radius of gyration R_g and the of bending stiffness was found. Unlike the papers of Cao et al. (2010), Ou et al. (2012), Singh et al. (2015), and Singh et al. (2016) made use of explicit water.

MD studies of insertion forces, inclusion free energy and osmotic pressure in coarse-grained polymer brushes were performed for particles of different shapes and sizes by de Beer et al. (2016) and different sizes, solvent conditions and degree of polydispersity by Merlitz et al. (2012) using the LJ and FENE potential and implicit solvent through the Langevin equation.

Brush configurations were studied by Li et al. (2018) who investigated the effect of counterion valence on the configuration of a spherical brush confined between two planes. The structure was probed by the end monomer distribution, ρ_m, brush thickness, R_g, and the radial distribution function g_np between brush monomers and counterions. Lastly, the mean square displacement of the counterions was found in directions parallel and perpendicular to the confining planes. The LJ, FENE and Coulomb potentials were applied.

Zhang and Xiang (2013) studied diffusion of a single free nanoparticle in a polymer brush through the particle-wall distance z_np, ρ_m, the lateral diffusion coefficient D_∥ and the force on the particle F_np. They found a decrease of D_∥ and increase of z_np with increasing grafting density. The simulation was set up using the Langevin equation and the LJ and FENE potentials. Csajka and Seidel (2000) used the same interactions, but also included the Coulomb potential in order to study polyelectrolyte brushes using MD. Instead of one diffusing particle, they had several counterions inside the brush. Their analysis included the counterion diffusion constants parallel and perpendicular to the plane, D_xy and D_z, which were both found to increase with the grafting density due to a decrease in the number of condensed counterions for denser brushes. While the work of Csajka and Seidel provides broad insight into polyelectrolyte brushes, it does not include a potential accounting for bending stiffness. Due to the semi-flexible nature of the constituents of the perineuronal nets (Richter et al., 2018), additional simulations are therefore required in order to understand diffusion in the PNNs. Furthermore, neither (Zhang and Xiang, 2013) nor (Csajka and Seidel, 2000) attempted to compare the diffusion constant to existing models of diffusion. Important insights might be gained from such models. Finally, estimating the effective resistance due to ion diffusion through a model polymer brush system may provide insights into the effects of PNNs on nerve signal propagation.

In this paper, we study diffusion in a charged, coarse-grained planar polymer brush by use of molecular dynamics as a first attempt at a model for diffusion in PNNs. The root mean square displacement of a single neutral or a single charged particle obtained through 1,000 system realizations is utilized to obtain the diffusion coefficient in directions perpendicular and parallel to the explicit substrate plane. Lennard-Jones, harmonic stretching and harmonic bending potentials are included in order to account for steric interactions, the stiffness of the bond and chain flexibility. Debye interactions account for screened electrostatic forces due water polarizability and dissolved ions in the liquid. The diffusion constants are found as a function of polymer chain spacing and compared to theory for a neutral particle. The diffusion constant for particles of various charges are then obtained and used to estimate the resistance of a 500 nm thick brush.

2. Model

Perineuronal nets consist of hyaluronan, chondroitin sulfate proteoglycans (CSPGs), Tenascin-R (Tn-R) and HAPLNs, the first acting as a scaffold and the latter two being cross-linkers (Richter et al., 2018). Hyaluronan (HA) and chondroitin sulfate (CS) are glycosaminoglycans (GAGs), which are unbranched, highly negatively charged sugar molecules of significant length (Richter et al., 2018). Although detailed information on the PNN structure is not known, it is hypothesized to look like Figure 1 (Richter et al., 2018). The nets form a cross-linked brush that is locally planar due to the relative sizes of the nets and the cell. As visible from Figure 1, several length scales are present in the system, namely the length of the chondroitin sulfate chains on the CSPGs, the length of the protein backbone of the CSPGs, the distance between chondroitin sulfate chains on the CSPGs, the distance between CSPGs on the hyaluronan chains, the length of the hyaluronan and the distance between neighboring hyaluronan chains. While a detailed model would be ideal, few of these PNN parameters are known. Instead of implementing a complicated, large-scale structure for which the exact structure is not known, the system was simplified to a planar glycosaminoglycan brush without attached CSPGs and cross-links in order to study diffusion in a GAG-rich structure. Investigating the effect of cross-linking and side-chains is left to future works.

FIGURE 1

Figure 1. A schematic of the perineuronal net structure.

As glycosaminoglycans are made up of repeating disaccharide units (Richter et al., 2018), larger systems can be modeled by coarse-graining the disaccharides into one bead instead of having explicit particles for every atom. One bead has length 1 nm, corresponding to a disaccharide unit of HA or CS, and the charge is −e just like in HA or unsulfated chondroitin. The disaccharide beads are linked together to form chains of N = 101 beads each, tethered to the membrane through two immobile beads. This length is short for membrane-attached HA, but a compromise was made in order to study sufficiently large system sizes. The chains proved sufficiently long to observe the most important effects of PNNs. The immobile beads form a 3 × 3 quadratic grid with lattice spacing d. To adhere to the GAG structure, chains are unbranched. Chains are kept flexible and of appropriate length through adjustable bond stretching and bond bending interactions, while the Lennard-Jones interaction prevents overlaps. The membrane is modeled by a rectangular lattice of immobile beads, whose LJ potential prevents other beads from passing through. A single unbonded particle is inserted close to the wall to probe diffusion in the brush. Periodic boundary conditions are applied in directions parallel to the substrate. A snapshot of the system is provided in Figure 2. The snapshot is rendered in OVITO (Stukowski, 2009), and two periodic images are included.

FIGURE 2

Figure 2. Snapshot of the modeled system for d = 15σ. The snapshot is rendered in OVITO with two periodic images in the x- and y- directions. Magenta beads-wall; coral beads-the two first beads of the chain (kept in place to ensure tethering to wall); blue beads-chain elements; yellow bead-diffusing particle (slightly enlarged for visibility). There is more than one yellow bead visible due to the periodic images, while only one particle is diffusing in the simulations.

3. Methods

3.1. Simulation setup

In order to capture properties of diffusion in the PNNs, molecular dynamics simulations of one particle in a system of tethered negatively charged polymer chains on a substrate were performed. The LAMMPS molecular dynamics package (Thompson et al., 2022) was used to run the simulations. Simulations were performed for 1,000 realizations of the system in order to gain appropriate statistics. Each realization had a different starting configuration of the particle and the chain beads.

The potentials used in the molecular dynamics simulations were selected to capture key properties of the PNNs. The Lennard-Jones potential was used to account for the spatial extent of the molecules. We have chosen the truncated 12–6 Lennard-Jones potential as it is common in the literature (see e.g., Merlitz et al., 2012). The potential $U_{\mathrm{\text{LJ}}}=\sum _{i<j}{u}_{\mathrm{\text{LJ}}}(r_{ij})$ works between all pairs of beads, including wall beads

$\begin{array}{ll}{u}_{\text{LJ}}(r)=4{ϵ}_{\text{LJ}}\left[{\left(\frac{{\sigma }_{\text{LJ}}}{r}\right)}^{12}-{\left(\frac{{\sigma }_{\text{LJ}}}{r}\right)}^6\right],r<r_{c,\text{LJ}},& (1)\end{array}$

where ϵ_LJ is the interaction strength and r is the distance between the beads. σ_LJ is the Lennard-Jones distance, which is related to the separation distance r_min of minimal potential energy by $r_{\text{min}}=2^{1/6}{\sigma }_{\mathrm{\text{LJ}}}$. r_{c, LJ} is the cutoff for the Lennard-Jones force.

As the repeating units of the GAGs have an approximately fixed size, we introduce a stretching potential following e.g., Zhang and Xiang (2013) and Csajka and Seidel (2000). Preliminary simulations showed that similar bond lengths were obtained using both the FENE potential (Kremer and Grest, 1990) and a harmonic stretching potential. We therefore chose a harmonic stretching potential as for example done by Horkay et al. (2020) and Kinjo et al. (2018). The harmonic stretching potential $U_{\mathrm{\text{bond}}}=\sum _{i<j}{u}_{\mathrm{\text{bond}}}(r_{ij})$ acts between every bonded pair of beads and is given as

$\begin{array}{ll}{u}_{\text{bond}}(r)=K_{\text{bond}}{(r-r_0)}^2,& (2)\end{array}$

where K_bond is the interaction strength, r₀ is the equilibrium bond length, and r is the distance between the bonded neighbors.

As the glycosaminoglycans comprising the majority of the PNNs are semi-flexible with a well defined persistence length (Richter et al., 2018), a harmonic bond bending potential was introduced as in Hehmeyer and Stevens (2005). The potential $U_{\mathrm{\text{bend}}}=\sum _{i<j<k}{u}_{\mathrm{\text{bend}}}({\theta }_{ijk})$ acts between every subsequent pair, ij, jk, of neighboring bonds

$\begin{array}{ll}{u}_{\text{bend}}(\theta )=K_{\text{bend}}{(\theta -{\theta }_0)}^2,& (3)\end{array}$

where K_bend is the interaction strength, θ₀ is the equilibrium angle, and θ is the angle between two subsequent bonds.

As the PNNs possess a highly negative charge, electrostatic interactions also need to be included. Furthermore, the cerebrospinal fluid contains dissolved ions that will screen charges in the extracellular matrix. The electrostatic interaction was therefore represented by the Debye potential. The Debye potential $U_{\mathrm{\text{Debye}}}=\sum _{i<j}{u}_{\mathrm{\text{Debye}}}(r_{ij})$ acts between all pairs of charged beads

$\begin{array}{ll}{u}_{\text{Debye}}(r)=\frac{C{q}_i{q}_j}{ϵr}{e}^{-\kappa r},r<r_{c,\text{Debye}},& (4)\end{array}$

where q_i and q_j are the bead charges, κ is the inverse of the Debye length (Nelson, 2014) and r_{c, Debye} is the cutoff for the interaction. C is an energy-conversion constant and ϵ is the dielectric permittivity. We assume that the Debye potential approximation is reasonable for modeling diffusion in interstitial fluids in the brain, with a Debye constant of approximately 1 nm (Syková and Nicholson, 2008).

The potential energy of the system is given by the sum of the stretching potential, the bending potential, the Lennard-Jones potential and the Debye potential

$\begin{array}{ll}U=U_{\text{bond}}+U_{\text{bend}}+U_{\text{LJ}}+U_{\text{Debye}}.& (5)\end{array}$

As incorporating explicit water molecules and ions is computationally expensive, the solvent was modeled implicitly. The Langevin equation is utilized as it mimics interactions with water by adding terms representing random collisions and drag effects of the fluid:

$\begin{array}{ll}m\frac{d\text{v}}{dt}=\text{F}(t)-\gamma \text{v}+\text{B}(t),& (6)\end{array}$

where B(t) is a random term (Lemons and Gythiel, 1997). The damping coefficient γ is defined as γ = 6πaη where a is the radius of the particle and η is the liquid viscosity.

Simulations were run for several types of systems. For the dynamic brush system (abbreviated as Dyn.), both the chains and the diffusing particle were considered dynamic particles and their time dynamics were calculated using molecular dynamics. For the static brush system (Stat.), the chains were modeled as static configurations by first equilibrating the chains, then fixing the chains while allowing the free particle to diffuse. This allowed us to compare the static and the dynamic system to provide insight into the role of geometry and flexibility on the diffusional properties of the system. In addition, simulations were performed on a dynamic brush system without bending terms (No stiff.) and in a static system with completely straight chains (Straight). Interactions between all particles, static or dynamic, were included in force and energy calculations.

Two types of grid geometries of the chain tethering points were studied: a quadratic and a hexagonal grid. However, preliminary simulations showed no systematic difference in the measured properties for the two grid types as illustrated in Supplementary Figure S1. We therefore limited the study to quadratic grids only.

A 3 × 3 quadratic grid was deemed sufficient as periodic boundary conditions were applied and no long-range forces were present. To verify that this assumption holds, the diffusion constants were found for systems of size 3 x 3, 4 x 4, 5 x 5, 6 x 6, and 9 x 9 for spacing d = 3σ. The diffusion constants agreed within the standard deviation. Furthermore, the sheer number of realizations for each grid spacing d should ensure that spatial heterogeneities are taken into account.

3.1.1. Parameter selection

The persistence length of a polymer characterizes its bending stiffness and is defined as the separation distance for which the correlation of bond vectors have fallen by 1/e. Its mathematical expression may be given as follows Kamerlin (2017):

$\begin{array}{ll}\langle cos\theta (s)\rangle =e^{-s/l_p},& (7)\end{array}$

where θ(s) is the angle between two bonds separated by the distance s along the polymer and l_p is the persistence length. For distances smaller than l_p, the bond orientations display a significant correlation (Kamerlin, 2017).

The parameters of the polymer model were tuned to yield a persistence length of 10 nm in accordance with Bathe et al. (2005), resulting in interaction parameters of K_bond = 140.2 and K_bend = 14.02. Lennard-Jones units were used, and the unit length σ was set to the length of a disaccharide unit, that is the bond length, so that σ = r₀ = 1 nm (Richter et al., 2018). The mass of the disaccharide unit is set to unity. The equilibrium angle was set to θ₀ = π in order to recreate the relatively linear nature of glycosaminoglycans.

In the Lennard-Jones interaction, the cutoff r_{c, LJ} is set to 1.122σ_LJ to coincide with the energy minimum at $r_{\mathrm{\text{min}}}=2^{1/6}{\sigma }_{\mathrm{\text{LJ}}}$, corresponding to good solvent conditions. The Lennard-Jones distance σ_LJ should be sufficiently large to avoid overlap between monomers. This is accomplished by ensuring that σ_LJ is equal to the length of a disaccharide unit, which means that σ_LJ = σ for chain beads. In the simulations the diffusing particle is the same size as the chain beads, σ_LJ = σ, except when finding D_bulk as a function of the particle radius a = σ_LJ/2. The interaction parameter ϵ_LJ = 0.73 in LJ units is derived from the carbon-carbon interaction parameter of 0.15 kj/mol that (Bathe et al., 2005) used in their coarse-graining of GAGs.

The chain beads hold a charge of −e each, corresponding to the charge of the repeating disaccharide unit of HA and unsulfated chondroitin (Richter et al., 2018). The substrate is neutral while the diffusing particle may be charged or neutral. The Debye length was set to match that of the brain extracellular fluid, namely 1/κ = σ = 1 nm (Syková and Nicholson, 2008), and the dielectric permittivity was set to match that of water. A cutoff of $r_{c,\mathrm{\text{Debye}}}=3{\kappa }^{-1}=3\sigma$ was implemented in agreement with Bathe et al. (2005).

4. Results

4.1. Anisotropic diffusion

A well-known result in physics is that the diffusion constant D_R in an isotropic system is related to the mean square distance 〈ΔR²(t)〉 traveled by the diffusing particles. In three dimensions, it takes the mathematical form of

$\begin{array}{ll}\langle \Delta R^2(t)\rangle =\langle \Delta x^2(t)+\Delta y^2(t)+\Delta z^2(t)\rangle =6D_{R}t,& (8)\end{array}$

where Δx²(t), Δy²(t) and Δz²(t) are the squared displacements in the x-, y- and z directions at time t, and ΔR²(t) is the total displacement.

The presence of a wall limits transverse motion, resulting in an anisotropic diffusion constant. Anisotropic diffusion has also been observed in brain areas with high concentration of PNNs (Morawski et al., 2015). In order accurately characterize diffusion in such systems, the mean square distance needs to be separated into contributions in directions parallel and perpendicular to the wall, as done by Carbajal-Tinoco et al. (2007)

$\begin{array}{ll}\langle \Delta x^2(t)+\Delta y^2(t)\rangle =4D_{\parallel }t& (9)\end{array}$

$\begin{array}{ll}\langle \Delta z^2(t)\rangle =2D_{\perp }t.& (10)\end{array}$

Here, D_∥ and D_⊥ are the diffusion constants in directions parallel (xy-plane) and perpendicular (z-axis) to the wall, respectively. Note that Equations (9) and (10) are valid in the diffusive regime, i.e., only for time scales longer than the initial ballistic regime. Strictly speaking, the diffusion constant of the particle should be a function of the particle-wall distance due to the anisotropy posed by the substrate. However, since the overall diffusive properties inside the nets are of interest, this is considered a detail outside the scope of this paper.

To ensure that the diffusion constant is characteristic for the brush, two versions of the mean square displacement (MSD) curve were found. The first curve is the MSD for the particle at all times regardless of whether or not it is still in the brush. At each point in time, the second curve is the MSD of particles still inside the brush. We disregard the initial part of the graphs, where inertial effects are prominent and the MSD curves are non-linear. The fits to Equations (9) and (10) were performed on an interval where the two MSD curves had similar slopes. In order to obtain error estimates, the MSD curves were divided into groups of 100 so that ten estimates of D_∥ and D_⊥ could be obtained using NUMPY's polyfit function (Harris et al., 2020). The average and standard deviation were calculated from these ten values. A plot illustrating the procedure is given in Supplementary Figure S2.

The free particle is subject to interactions with chain- and substrate beads in the system. For interacting particles, two diffusive regimes occur. On shorter time scales, the particle moves close to a local energy minimum in the energy landscape (Dhont, 1996). As time progresses, the particle passes many energy barriers and visits many local minima. Additionally, changes in brush configuration alter the potential energy landscape, creating new minima for the particle to visit. The short-time and long-time diffusion constants are therefore different. Between these time scales, where only a few minima has been visited, the MSD is non-linear (Dhont, 1996). Our MSD plots did not display such a non-linear interval while the particle was inside the brush, indicating that long-time diffusion was not reached. Note that particles exiting the brush puts an upper limit on the time scale of diffusion.

Figure 3 shows the bulk diffusion constant D_bulk as a function of the radius a of the diffusing particle found by Equations (8)–(10). The bulk diffusion constant was found by tracking the movement of a single particle subject to only Equation (6), where γ depends on a. The system is otherwise empty, meaning that F=0 in Equation (6) in this instance. We see a sharp decrease of D_bulk with increasing a. The diffusion constants in all directions are within the error margins of each other, as expected due to the isotropic nature of the bulk system. The insert in Figure 3 shows a log-log plot of D_bulk vs. a. The curves appear linear except for some small deviations at the lowest value of a = 0.125 nm. The data is in agreement with the Stokes-Einstein equation D = kT/(6πηa) (indicated by the dashed line) within the standard deviation (Edward, 1970).

FIGURE 3

Figure 3. Diffusion constant D vs. radii a of a free particle in bulk. D_R - total diffusion constant; D_⊥-diffusion constant in the z-direction; D_∥-diffusion constant in the xy-plane; kT/γ-the theoretical value of D. Standard deviations are indicated by shaded regions. Insert: A log-log plot of the same data.

4.2. Brush configuration vs. spacing d

Figure 4 shows system configurations for spacings of d = 3σ, d = 10σ and d = 50σ. For the smallest values of d, such as d = 3σ in Figure 4A, the chains are far straighter than for larger spacings such as d = 10σ and d = 50σ. Table 1 shows the potential energy of the neutral diffusing particle and reveals that at small d, the particle experiences both stronger and more frequent interactions with other beads in the system. Also included in Table 1 is the average potential energy of the chain beads. The chain beads experience similarly increased Debye and Lennard-Jones interactions at small d, forcing them into straighter configurations in order to lower the energy. For d < 2σ, visualization of the trajectories revealed an initial period of random motion until the diffusing particle starts moving back and forth within a small volume in the brush. A video of a trajectory for d = 1.5σ is provided in Supplementary Video S1. Keeping in mind that the Lennard-Jones interaction has a minimum at $2^{1/6}{\sigma }_{\mathrm{\text{LJ}}}$ and that σ_LJ = σ for both particle and beads in these simulations, the particle will experience significant forces from multiple chains propelling it forwards until it happens to encounter a cavity within the brush. If the particle velocity is sufficiently small or the cavity is sufficiently large, the particle will not be able to move forward, but will be trapped inside the cavity. At longer time scales, the particle would typically have been able to escape this volume, displaying a reduced but non-zero diffusion constant in the long-time diffusion regime (Cai et al., 2011). Such a diffusion regime was not reached for our simulations.

FIGURE 4

Figure 4. Brush configurations for different spacings d. (A) d = 3σ, (B) d = 10σ and (C) d = 50σ. The chains are much straighter in A than in C. Magenta beads-wall; blue and coral beads-chain beads; yellow bead-free particle (slightly enlarged for visibility). The figure is rendered in OVITO.

TABLE 1

Table 1. Potential energy of the free particle and the brush.

As d increases and beads on different chains get farther apart, the magnitude of the interactions between beads on different chains decreases, as seen in Table 1, and the chains take on a more curled configuration, as seen for d = 10σ in Figure 4B. As d is increased further, interchain interactions will rarely occur. The system is then said to be in the mushroom configuration, which is typically distinguished from the brush configuration (Kim et al., 2015).

In summary, the system is shown to enter three regimes: The extremely dense brush for d < 2σ, in which neutral particles get trapped at our time scales, the mushroom regime for large d, in which chains are far apart and we expect diffusion to be similar to bulk diffusion, and the brush regime for intermediate values of d. Simulations for large values of d have been performed in order to verify that the system tends to bulk behavior in the limit of large d. Intermediate d-values are of primary interest, as they yield systems similar to that of perineuronal nets. We therefore select d = 2σ as the lower limit for d, ensuring that we study a range of d-values for which the particle is sufficiently mobile in most of the simulations. The aptness of this limit can be seen by comparing the mean square displacement in Supplementary Figures S3, S4.

Tethered polymers are in the brush configuration for d < 2R_g (Kim et al., 2015), where R_g is the radius of gyration, which for a polymer chain is defined as Kamerlin (2017):

$\begin{array}{ll}{R}_g=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N}{({\text{r}}_i-\bar{\text{r}})}^2},& (11)\end{array}$

where N is the number of monomers in a chain, r_i is the position of the ith monomer and $\bar{\text{r}}$ is the center of mass. Figure 5 shows 〈R_g〉 vs. d for the dynamic brush plotted in the form of a phase diagram, where the brush and mushroom regimes are indicated by the background color. We found 〈R_g〉 by averaging R_g in Equation (11) over all chains and realizations for the last time frame of each simulation. Based on Figure 5 we find that the system is in the brush configuration for d ∈ [2σ, 25σ], whereas for larger values of d the system is in the mushroom configuration. Here, we have focused on the brush regime which is most similar to PNN structure, and have primarily analyzed systems with d < 25σ.

FIGURE 5

Figure 5. R_g vs. d for the dynamic brush. The grafted polymers are in the brush configuration for d < 2R_g and in the mushroom configuration for d > 2R_g (Kim et al., 2015), as indicated by the colored background. The standard deviation is indicated by the shaded blue area.

To compare the brush configuration to theoretical predictions, the average height of the brush L_z was found for each d by extracting the maximal z-coordinate of the chain beads at the last time frame and averaging over all realizations. The brush height as calculated by mean field theory scales with brush spacing as $L_z=A{d}^{-2/3}$ (Attili et al., 2012), where A is a constant that does not depend on d. The data points were fit to theory using the Levenberg-Marquardt algorithm as implemented in SCIPY's optimize.curve_fit function (Virtanen et al., 2020). The height of the brush agreed fairly well with the theory, as shown in Supplementary Figure S5 and Supplementary Table S1, where fits to theory are included as dashed lines. The height of the brush without bending potential energy showed the best agreement with the theory, while the systems with bending stiffness also displayed an acceptable agreement. The system types are as defined in Section 3.1.

4.3. Diffusion constant D vs. grid spacing d

4.3.1. Dynamic brush

Figure 6A shows D(d)/D_bulk for the dynamic and the static brush. This ratio converges to unity when the lattice spacing d increases, which is as expected since the system effectively approaches bulk as d becomes large. The decreasing values of D as d becomes smaller are due to brush interactions, which become more dominant as the density of the brush increases. At d < 2σ frequent entrapment yields an effective diffusion constant which is close to zero.

FIGURE 6

Figure 6. D/D_bulk vs. d in different directions for a neutral particle in dynamic and static brushes. (A) D/D_bulk for d ∈ [1σ, 100σ]. D/D_bulk ≈ 0 for d < 2σ, then increases with d and eventually converges to 1. (B) D/D_bulk for d ≤ 10σ. Diffusion constants for straight, immobile chains and chains without a bending term are included. Blue line, D_⊥; green line, D_∥. Bold colors, dynamic brush; light colors-static brush. Red, D_⊥ for the straight system; light red-D_⊥ for the system without bending stiffness. Black, D_∥ for the straight system; light gray, D_∥ for the system without bending stiffness. Standard deviations are indicated by shaded regions.

Figure 6B shows the data in Figure 6A for up to d = 10σ, with the additional systems of totally straight chains and chains without a bending term included in the plot. The curves indicate that D_∥ < D_⊥ from d = 2σ and up to about d = 10σ. As seen by Figure 4A, the chains fall into straighter conformations when the spacing is small. This is, as can be seen from Table 1, due to increased forces between the beads composing the chains. The pore space between the chains then take on a shape similar to channels oriented in the z-direction, opening for more rapid motion in the perpendicular direction. However, the proximity of the brush will still lead to collisions, lowering D_⊥ compared to bulk. For motion parallel to the substrate, there are no such open channels, leading to a larger reduction of D_∥ compared to D_⊥. In summary, we expect D_∥ < D_⊥ for smaller values of d due to the anisotropic geometry of the system. Looking at Figure 4B, the chains appear to generate a more disordered structure for d = 10σ, and do not display similarly visible channels as they do for d = 3σ in Figure 4A. Inspecting Figure 6B further, we see that the difference between D_⊥ and D_∥ is larger for small d, where the channel-like formations are more prominent.

4.3.2. Static brush

Figure 6 also shows D(d)/D_bulk for the static brush. As for the dynamic brush, D_∥ < D_⊥ for small values of d, and both D_⊥ and D_∥ converge to unity as d → ∞. Figure 6B shows that D_{⊥, stat.} < D_{⊥, dyn.} for smaller d, although D_⊥ of the two system types agree within the standard deviation. This miniscule increase in the diffusion constant may be due to the flexibility of the dynamics chains, allowing the chains to sway due to the motion of the diffusing particle. This small difference in the diffusion constants of the static and dynamic brushes indicate that the dynamic flexibility of the chains do not play a major role in governing diffusion of free particles within the brush. The systems have similar geometry, as the chains have been made immobile only after equilibration. We therefore conclude that the geometry and not the dynamics of the brushes determines the diffusion coefficients.

4.3.3. Other system types

Also included in Figure 6B are the diffusion constants of a free particle in a system where the bending potential energy (as defined by Equation (3)) is removed, and a system where the chains are straight and immobile. Due to the chain straightness in the latter system, we expect the channels to be more pronounced and that motion in the perpendicular direction will be less restricted. We therefore expect D_{⊥, straight} to be larger than D_⊥ for any of the other systems, as demonstrated by Figure 6B. For simulations where the bending stiffness term is removed, the chains curl up more, yielding a slightly lower D_⊥, as seen from the coral curve (marked by X's) in Figure 6B.

In the z-direction the geometry varies from clearly defined channels for the straight chains to mildly obstructed channels for the static and dynamic chains. In the direction parallel to the wall, along the xy-plane, there are few clear channels, but rather a forest of obstacles. Like the dynamic system, neither the static system nor the system without bending interactions exhibit ordered obstacles in the xy-plane. Although the straight system is more ordered than the other systems, it still exhibits the forest-like geometry in the xy-plane that all the system types display. We therefore expect that every system type possess similar values of D_∥/D_bulk. Figure 6B reveals that the D_∥/D_bulk curves of the various system types all agree well within their standard deviations, unlike the curves of D_⊥/D_bulk.

4.4. Effective diffusion behavior

4.4.1. Models for effective diffusion behavior

Many theoretical and experimental efforts strive to connect the diffusion constant to the porosity ϕ (Shen and Chen, 2007), that is, the fraction of space not occupied by solid material. In the simulated systems, the chain beads constitute the solid material, while the porosity is effectively the solvent volume fraction.

The porosity of our system can most easily be estimated through the total volume of the beads, V_beads, and the volume of the system, V_system. Using the volume of a sphere, $V_{\mathrm{\text{beads}}}=N_{\mathrm{\text{chains}}}{N}_{\mathrm{\text{bead}}}·\frac{4}{3}\pi a_{\mathrm{\text{f}}}^3$, where N_chains is the number of chains, N_bead is the number of beads per chain and a_f is the radius of the beads. $V_{\mathrm{\text{system}}}=N_{\mathrm{\text{chains}}}{d}^2{L}_z$ was then found, where L_z is the average height of the chains. The porosity ϕ of the brush is hence

$\begin{array}{ll}\varphi =1-\frac{\frac{4}{3}\pi a_{\text{f}}^3{N}_{\text{bead}}}{d^2{L}_z},& (12)\end{array}$

However, since there are many porous systems with the same porosity, but with very different geometries, we expect the diffusion constant to depend on further properties of the porous system geometry. In our simulations, this is evident from comparing D_⊥/D_bulk of the different system types, as seen in Supplementary Figure S6, where the system of straight chains displays a higher diffusion constant than the others for a wide range of ϕ's.

One property that affects the diffusion is the tortuosity τ, a measure of the degree to which diffusing particles are prevented to move in straight lines due to obstacles posed by the porous media. We may define τ as the path length Δl traversed by the diffusing particle divided by the corresponding linear distance Δx (Shen and Chen, 2007).

$\begin{array}{ll}\tau =\frac{\Delta l}{\Delta x}.& (13)\end{array}$

Using this definition, τ yields a correction to the bulk diffusion constant D_bulk that reads

$\begin{array}{ll}{D}_{\text{eff}}=\frac{D_{\text{bulk}}}{{\tau }^2},& (14)\end{array}$

where D_eff is the effective diffusion constant in the porous medium. Shen and Chen (2007) list τ² as a function of ϕ for several previous studies.

Having obtained the porosity ϕ, the diffusion constant D(d) can be compared to theory in order to characterize diffusion. Due to the anisotropic nature of the system, D_⊥ and D_∥ might be best described by different diffusion models. These models can subsequently be used to offer an approximation of D_⊥ and D_∥ for values of d that have not been simulated, but are within the range of interest. Furthermore, we can generalize from one bead size to another by using the parameters obtained in the fit in combination with Equation (12). Ultimately, this might prove useful in obtaining other net quantities from D, such as the capacitance.

The hyperbola of revolution listed in Shen and Chen is a system whose pore space comprises several consecutive hyperbolas (Petersen, 1958). Its tortuosity τ is defined through

$\begin{array}{ll}{\tau }^2=2-\varphi .& (15)\end{array}$

Another expression for τ reads

$\begin{array}{ll}{\tau }^2=(3-\varphi )/2.& (16)\end{array}$

This relation holds for various systems, such as for a dilute system of spheres (Akanni et al., 1987), spheres of various sizes that are allowed to overlap (Akanni et al., 1987) and randomly packed spheres of different radii (Neale and Nader, 1973). These systems are dubbed “ordered packings” by Shen and Chen (2007).

A tortuosity model frequently applied to heterogeneous catalysts was developed by Beeckman (1990), Yang et al. (2014), and Ghanbarian et al. (2013), where

$\begin{array}{ll}{\tau }^2=\varphi /\left[1-{(1-\varphi )}^{1/3}\right].& (17)\end{array}$

The cation-exchange resin membrane consists of cross-linked polymers with uniformly distributed negative charges (Mackie and Meares, 1955). Its tortuosity was derived by Mackie and Meares (1955) using a random walk lattice model. A site in the model can either be occupied by polymer, solution or the walker. The walker is forced to move around polymer sites (Mackie and Meares, 1955), which it will encounter with a probability of p_first = ϕ_p = 1 − ϕ, where ϕ_p is the polymer fraction. Upon encountering a polymer site, there is a set probability p_next that the next site in the walk will also contain a polymer. This probability is chosen to reflect the properties of space in the cross-linked polymer and is found to be $p_{\mathrm{\text{next}}}=\frac{1}{2}(2-\varphi )$ by Mackie and Meares (1955). The probability that the walker will have to move two extra steps to get past the polymer is then $p_1=(1-\varphi )\left[\frac{1}{2}(2-\varphi )\right]$, the probability of three extra steps is $p_2=(1-\varphi ){\left[\frac{1}{2}(2-\varphi )\right]}^2$ and so on. The total length Δl the walker will have to move is related to the path length in bulk Δx through Mackie and Meares (1955):

$\begin{array}{ll}\tau =\frac{\Delta l}{\Delta x}=1+(1-\varphi ){\displaystyle \sum _{n=0}^{\infty }}{\left[\frac{1}{2}(2-\varphi )\right]}^n,& (18)\end{array}$

where n is an index running over all the possible times the walker can hit polymer elements consecutively. Per definition, the pore fraction ϕ cannot exceed one, and for ϕ = 0 the system would only consist of polymer. Hence 0 < ϕ < 1, which means that $\frac{1}{2}(2-\varphi )<1$ and that the geometric sum in Equation (18) converges. Performing the sum and adding the first term in Equation (18) yield (Mackie and Meares, 1955):

$\begin{array}{ll}\tau =\frac{2-\varphi }{\varphi }.& (19)\end{array}$

Combining Equations (14) and (12) with Equations (15)–(17) and (19) allows for a comparison between data and theories.

4.4.2. Modified model for effective diffusion behavior

The probability p_next in Mackie and Meares was tailored to the cation-exchange resin membrane. However, the probability p_next can be modified to better capture the geometry of other systems such as the PNN system we study.

Keeping the probabilities p_first and p_next general for the sake of brevity, the probability that the walker will have to move two extra steps to get past the polymer is p₁ = p_firstp_next, the probability of three extra steps is $p_2=p_{\mathrm{\text{first}}}{p}_{\mathrm{\text{next}}}^2$ and so on. We thereby arrive at a more generalized option for Equation (18).

$\begin{array}{ll}\tau =1+p_{\text{first}}{\displaystyle \sum _{n=0}^{\infty }}{p}_{\text{next}}^n,& (20)\end{array}$

where n is an index running over all the possible times the walker can hit the polymer consecutively. As p_next is a probability and p_next = 0 corresponds to bulk behavior, 0 < p_next < 1. Performing a geometric sum on the last term in Equation (20) yields:

$\begin{array}{ll}\tau =1+\frac{p_{\text{first}}}{1-p_{\text{next}}},& (21)\end{array}$

which offers a generalization to Equation (19).

Several different probabilities p_next were tested against the diffusion coefficients the polymer brush, and the best p_next was selected as a custom model. It was derived by assuming that if the walker hits a polymer section, it has a probability of p_next = k + (k − 1)(1 − ϕ) of not being able to make the next move because a neighboring site is occupied. Here, k is a tunable parameter. Since the solid is solely made up of polymer, the polymer fraction 1 − ϕ should equal the probability to first hit an obstacle, p_first. The relation p_first = 1 − ϕ from Mackie and Meares (1955) should therefore hold for all systems and is kept as is. Insertion of p_first and p_next into Equation (21) yields an expression of the tortuosity:

$\begin{array}{ll}{\tau }^2={\left[1+\frac{1-\varphi }{1-k-(k-1)(1-\varphi )}\right]}^2.& (22)\end{array}$

This model will be referred to as the custom model. Combining Equations (14) and (12) with Equation (22) allows for a comparison between data and model.

4.4.3. Power law model

Judging by the shape of the curves in Figure 6, we may approximate D/D_bulk as a power law:

$\begin{array}{ll}\frac{D}{D_{\text{bulk}}}=1-d^{-n}c,& (23)\end{array}$

where n and c are tunable parameters for this model.

4.4.4. Comparison with data

As none of the models listed in Shen and Chen (2007) closely resemble the brush systems we study, the bead size a_f was left as a free parameter for the curve fit using the Levenberg-Marquardt algorithm through SCIPY's optimize.curve_fit function (Virtanen et al., 2020), as the theoretical models were not tailored to the brush and some discrepancy should be expected. The optimal parameter values for each fit to the diffusion constant of a neutral particle in the dynamic brush and the static brush from Figure 7 is given in Table 2, together with the relative differences between the parameters in the parallel and perpendicular directions. The parameter k in the custom model was found by the Trust Region Reflective algorithm as implemented in SCIPY and n and c for the power law was found by the Levenberg-Marquardt algorithm as implemented in SCIPY. The results are listed in Table 3 for all directions for the static and dynamic brush. The brush regime is of interest for the fits since this configuration should correspond to developed PNNs (Richter et al., 2018). Furthermore, the PNNs are reported to be quite dense (Deepa et al., 2006). The range d = 2 − 10σ was chosen to satisfy as closely as possible both these criteria. Figure 7 shows D_∥/D_bulk and D_⊥/D_bulk for the dynamic and static brushes together with fits. For readability purposes, only three models are selected for plotting: The hyperbola of revolution, which returns the same curves as the ordered packings, the custom model and the power law model. The complete set of fits including different custom models are shown in Supplementary Figures S7–S10 and specified in Supplementary Tables S2, S3. The cation-exchange resin model proved a poor fit in all instances, while the heterogeneous catalyst was slightly outperformed by the power law for D_⊥, while providing a poor fit for D_∥.

FIGURE 7

Figure 7. Diffusion constants as a function of d together with fits. Solid line-D_∥/D_bulk and D_⊥/D_bulk vs. d for dynamic and static brushes. (A) D_∥/D_bulk for the dynamic brush. (B) D_∥/D_bulk for the static brush. (C) D_⊥/D_bulk for the dynamic brush. (D) D_⊥/D_bulk for the static brush; D, data; hr, Hyperbola of revolution fit; c, Custom model fit; pl, Power law fit. Standard deviations are indicated by shaded regions. Parameters of the fits are listed in Table 2.

TABLE 2

Table 2. Parameters a_f of the fits to models listed in Shen and Chen (2007) for the dynamic and static brush.

TABLE 3

Table 3. Parameter k of the fit to the custom model and parameters n and c of the fit to the power law.

Figure 7A reveals that the power law of Equation (23) offers a poor fit to D_∥/D_bulk for the dynamic brush. A better fit to the data is provided by the hyperbola of revolution. The custom model overall offers the best fit to the data, although it starts to deviate slightly for larger values of d. This might be an artifact of the data set, since there are relatively fewer data points for large d compared to small d.

In Figure 7B D_∥(d)/D_bulk for the static brush is plotted together with fits to the different models. The best parameter fits for the static brush are given in Table 2 along with the relative differences between the results in the parallel and perpendicular directions. The parameters a_f are consistently larger for parallel motion, with the exception of the cation-exchange resin membrane model. The custom model and the power law provide fits of comparable quality, with both models deviating from the data set for larger d's and the custom model providing a better fit at smaller d's. Both fits mostly stay within the standard deviation. The hyperbola of revolution more closely aligns with the data set. Considering the unrealistically large value of a_f for the hyperbola of revolution in Table 2, the custom model proves the best model overall for D_∥(d)/D_bulk in the static brush.

As seen by comparing Figures 7A,C, D_⊥/D_bulk for the dynamic brush generally has larger uncertainties than D_∥/D_bulk for d ∈ [2σ, 10σ]. The fitted parameters are therefore more uncertain. Both the custom model and the hyperbola of revolution perform poorly on this system, overestimating the diffusion constant for larger d and slightly underestimating it for smaller d. The power law more closely aligns with the data points, staying within a standard deviation while retaining a similar curvature. Similar results are found for D_⊥/D_bulk for the static brush, as shown in Figure 7D.

It is worthwhile to note that the theoretical models were originally developed for describing diffusion in three-dimensional space. By comparing to two-dimensional diffusion in the case of D_∥ and one-dimensional diffusion in the case of D_⊥, the models lose some of their interpretive qualities. Instead, they provide a heuristic mathematical framework. The power law often takes on this role, as it is abundant in nature (Gheorghiu and Coppens, 2004).

4.5. Relation between D_⊥ and D_∥

Both diffusion constants, D_⊥ and D_∥, approach D_bulk when the spacing d becomes large. However, the way they approach D_bulk can provide insights into the role of anisotropy in the system. To address this, we plot the ratio D_⊥/D_∥ as a function of spacing d in Figure 8. For an isotropic system, we would expect the diffusion constants to depend on d, but the ratio to be constant. Figure 8 shows that D_⊥/D_∥ approaches unity for large d. The behavior of D_⊥/D_∥ is well approximated by a power law, $D_{\perp }/D_{\parallel }=A{d}^{-l}+1$, as illustrated by the plot of the fitted function in Figure 8 and in the log-log plot insert. We performed similar analysis for all systems (see Supplementary Figures S11–S13) and the results are shown in Table 4. The range of spacings, d, over which the systems display significant effects of anisotropy is within the range relevant for PNNs (Deepa et al., 2006; Richter et al., 2018). We therefore expect anisotropic effects in diffusion to be important for PNNs.

FIGURE 8

Figure 8. D_⊥/D_∥ vs. d for the dynamic brush, together with a fit to the power law Ad^−l + 1. Insert: A log-log plot of D_⊥/D_∥ − 1 together with Ae^−d/l. The fit is performed in the range d ∈ [2.5σ, 25σ]. Standard deviations are indicated by shaded regions. Insert: A log-log plot of D_⊥/D_∥ − 1 and Ad^−l. Note that D_⊥/D_∥ − 1 fell below zero for d = 15σ and d = 25σ so that these points are not visible in the log-log plot.

TABLE 4

Table 4. Parameters of the fits of D_⊥/D_∥ to f(d) = Ad^−l + 1 for the different system types.

4.6. Diffusion vs. particle size

Simulations were performed for a diffusing particle of various radii a in the dynamic brush. The resulting perpendicular diffusion constants were divided by D_⊥ for a = 0.5 nm to illustrate how diffusion is affected by particle size.

Scaling theories of diffusion in polymer liquids predict that the diffusion constant scales with particle diameter σ_LJ as D ∝ 1/σ_LJ when σ_LJ < ξ (Cai et al., 2011), where ξ is the correlation length. ξ is defined the average distance from one chain bead to the closest bead on another chain (Cai et al., 2011). We found ξ by looping over all beads for the last time frame of the simulation and finding the closest bead on another chain. The average was taken over all beads and realizations. A plot of ξ vs. d is presented in Supplementary Figure S14, revealing that σ_LJ < ξ for the majority of spacings studied. D_⊥/D_{⊥, a = 0 nm} is therefore compared to σ_{LJ,a = 0 nm}/σ_LJ = 1/2a.

Results showed an acceptable agreement with theory, as shown in Figure 9A. D_⊥ increases with decreasing a, as expected from Equation (6). The exception is for particle diameters close to d, i.e., when the brush is dense compared to particle size. Inspection of trajectories reveals occasional extrusion of the particle from the brush for such parameter combinations, probably due to unfavorable initial conditions for these large particles.

FIGURE 9

Figure 9. Diffusion properties for an uncharged particle of different radii a in a dynamic brush. (A) D_⊥/D_{⊥(a = 0.5nm)} vs. d, where D_{⊥(a = 0.5 nm)} is D_⊥ for a= 0.5 nm. The dashed lines indicate the expected scaling of the diffusion constant as presented in Cai et al. (2011). (B) Diffusion times through the brush vs. d for a brush height of h = 500 nm. The curves were found through application of Equation (10).

We applied (Equation 10) to find the corresponding diffusion time through a brush of height h = 500 nm as a function of spacing d, as shown in Figure 9B. The diffusion time stays within 0.4 μs for a ≤ 0.5 nm.

4.7. Diffusion of charged particles

4.7.1. Diffusion constants

Ionic transport and the diffusion of molecules with more complex surface charge distributions are important in PNNs. We therefore address the diffusion of charged particles, which may represent an ion or a charged molecule, through the brush as a model for transport through PNNs. Figure 10A shows D_∥/D_bulk vs. d for diffusing particles with different charges. We see that the diffusion constant is very small for positively charged particles compared to negatively charged and neutral particles due to electrostatic interactions. Positively charged particles are therefore subject to an intermittent and slow transport mechanism, jumping between trapping sites that are close to the chain beads (see Supplementary Videos S2–S4) in a process similar to cations binding to negatively charged sites on GAG chains in the PNNs (Morawski et al., 2015). Negatively charged diffusing particles experience repulsive electrostatic forces from the chain beads, which effectively increases the particle interaction size, and therefore restricts the available pore space and lowers of the diffusion constant.

FIGURE 10

Figure 10. D/D_bulk vs. d for a diffusing particle of charge q. (A) D_∥/D_bulk vs. d. (B) D_⊥/D_bulk vs. d. Standard deviations are indicated by shaded regions.

Figure 10B show that the behavior of D_⊥/D_bulk for charged particles is similar to that of D_∥/D_bulk. The increased value of D_⊥ for small d (d < 3.5σ) is likely an artifact from the simulation procedure: If the initial position of the charged particle had been outside of the brush, the repulsive electrostatic forces would have most likely have prevented it from entering the brush. Lower concentrations of anions compared to cations have been observed in GAG-rich environments (Morawski et al., 2015), so it seems reasonable that the brush should exclude anions at sufficiently small d. Overall, the diffusion constants are reduced by a factor of two when the charge was increased from q = 0 to q = −e.

4.7.2. Conductivity, resistivity, and resistance

The conductivity can be computed from the diffusion constant of ions in solution through

$\begin{array}{ll}{\sigma }_e=\frac{F^2}{RT}{\displaystyle \sum _x}{D}_x{z}_x^2\left[x\right],& (24)\end{array}$

where σ_e is the electrical conductivity, F is Faraday's constant, R is the gas constant, D_x is the diffusion constant of ionic species x, z_x is the valency of the ion and [x] is its concentration (Qian and Sejnowski, 1988). To find the conductivity σ_e, we used concentrations from Raimondo et al. (2015), adding together concentrations of ions of the same valency for use in Equation (24).

Figure 11A shows the electrical conductivity σ_e as a function of d across the brush. Here, we have excluded small values of d since these systems are unrealistic for charged particles as argued above.

FIGURE 11

Figure 11. Electrical properties of the brush as a function of d. (A) Electrical conductivity σ_e(d) through the brush. (B) Resistance R vs. d for different cell radii r_neuron assuming a brush height of h = 500 nm. Standard deviations are indicated by shaded regions.

The resistance R is connected to the resistivity ρ_e = 1/σ_e by

$\begin{array}{ll}R=\frac{{\rho }_{e}l}{A}& (25)\end{array}$

where l is the length of the conducting material and A is its cross-section. In cases where A varies along the length of the conductor, the equation should be integrated in order to yield correct results. PNNs wrap around the soma and proximal dendrites of a neuron (van 't Spijker and Kwok, 2017). Assuming that PNNs take the form of a spherical shell of thickness h, the resistance of the nets become:

$\begin{array}{ll}R=\frac{{\rho }_e}{4\pi }\left(\frac{h}{r(r+h)}\right)& (26)\end{array}$

where r is the radius of the neuron. The resistance is inversely proportional to the diffusion constant, as seen from Equation (24), meaning that a lowered diffusion constant leads to an increased resistance.

Figure 11B shows the resistance R vs. d for a brush with height h = 500 nm for different cell radii r_neuron (Quan et al., 2013). We see that R decreases with increasing d, which is expected because as the brush becomes sparser, the average interactions between brush and particles decrease.

We found the effective resistance of the brush to be R ≃ 1 − 62 Ω. Input resistances of neurons are usually much higher, on the order of magnitude of 100 MΩ (Jouhanneau et al., 2018). With their small resistance, the brushes do not appear to pose a significant barrier to current of negatively charged ions.

5. Discussion

We have used D_∥/D_bulk and D_⊥/D_bulk to characterize diffusion inside a brush as model for transport in PNNs. The values of D_⊥/D_bulk in Figure 6 show that a charged brush limits diffusion in the transverse direction as mentioned by van 't Spijker and Kwok (2017), an effect that increases with increased density. Low values of D_⊥ in dense brushes supports that glycosaminoglycans block out molecules that may be potentially harmful to the cell (Suttkus et al., 2014; Morawski et al., 2015). Lateral diffusion, as illustrated by D_∥/D_bulk in Figure 6, exhibit a similar behavior.

Negatively charged free particles have a lower diffusion constant than uncharged particles, highlightning the ability of the net to act as a barrier to ions, as has been proposed by others Morawski et al. (2015); van 't Spijker and Kwok (2017). The drastically lowered diffusion constants of positively charged particles as shown in Figure 10 are to be expected as the positive charges experience an attraction to the negatively charged nets by electrostatic forces. However, this may only be a transient effect until the negative charge has been neutralized by trapped positively charged ions.

Equation (22) represented a good heuristic model for diffusion in this system. It is based on the generalized version of the Mackie and Meares (1955) model that we introduced in Equation (21), which could prove a useful model for the diffusion constant in a wide range of systems. The parameter p_next can be estimated for a specific geometry, providing a method to predict τ for given geometries.

The viscosity η in Equation (6) was assumed to be constant and equal that of bulk. However, when the diameter of the diffusing particle is larger than the correlation length ξ, an effective viscosity arises in polymer liquids (Cai et al., 2011). This viscosity is increased compared to bulk (Cai et al., 2011). Effective viscosities have also arisen in polymer brushes (Doyle, 1998). In the present study, the particle diameter σ_LJ < ξ for all d, and the same holds for all data points presented in Figure 9, as seen by Supplementary Figure S14.

Representing the perineuronal net structure by a charged planar brush is a strong simplification. As seen by Figure 1, perineuronal nets consist of several components that vary in shape and size, which are bound and cross-linked to form a complicated structure for which the details are currently not known. An important aspect of PNNs is the highly negative charge density of the GAGs (Richter et al., 2018). Since we have included a high charge density in our model, the diffusional behavior should still capture key aspects of motion into and out of the perineuronal nets.

Relatively small values of d were used in this study compared to what should be realistic for PNNs. Real-life PNNs will have side CSPGs attached, akin to side chains. The CSPGs have a certain spatial extent, increasing the effective volume of each HA chain in the nets. A model consisting of linear chains correspondingly underestimates the volume fraction of polymer inside the brush. Lowering d will therefore yield a more accurate porosity, which is an important concept in diffusion in porous media. Given that accurate measurements of grafting densities and CSPG spacing on HA in PNNs are not currently available, this approach seemed favorable over creating a detailed, more realistic model with uncertain parameters. We simulated systems with a wide range of spacings, but emphasis was placed on the range d = 2 − 10 nm for which diffusion was significantly different than in bulk. This range was chosen as PNNs are hypothesized to restrict transport (Morawski et al., 2015; van 't Spijker and Kwok, 2017; Fawcett et al., 2019).

Estimates of persistence lengths of each GAG vary somewhat, with persistence lengths as short as 4 nm being measured for HA (Bano et al., 2016). The persistence length is governed by the relative size of the bending interaction parameter K_bend in Equation (2). This force field term was absent in the system type with no stiffness, meaning that it should have a persistence length close to the length of the bead, that is 1 nm. The stiffness of the chains did not affect diffusion noticeably, as seen by comparing the diffusion constant in systems with bending terms (dynamic) and without (no stiffness) in Figure 6B. Hence, the exact value of l_p does not pertain significantly to the study of diffusion constants in these brushes.

PV interneurons, around which PNNs frequently enwrap, are well-known to have low input resistance (Ferguson and Gao, 2018; Jouhanneau et al., 2018). If the net resistance was large or comparable to the membrane resistance, the total resistance of PV interneurons would be large compared to other neurons, or in the very least not significantly lower. The low resistance of our simplified PNN model is hence in agreement with experimental results, but may have corrections due to the simplifications made.

The PNNs are reported to protect neurons against ions and molecules of different sizes (van 't Spijker and Kwok, 2017), so particles smaller than the chain beads will sometimes be of interest. The Lennard-Jones potential (Equation 1) is included to account for steric interactions, and σ_LJ and hence r_{c, LJ} should therefore depend on particle size. A smaller particle will thus be less influenced by the Lennard-Jones potential as it is within the cutoff for a smaller amount of time. Furthermore, the mass should be adjusted according to size, yielding a smaller particle mass. Taking both radius and mass into account, the damping term and the random term in Equation (6) will be lower. The motion of the particle will therefore be faster, and the diffusion constant in brush and in bulk will be a bit higher. This is in agreement with Figures 3, 9.

The charged particle is sensitive to the Debye interaction, which does not depend on particle size, only separation distance. When decreasing the particle radius, the Lennard-Jones interaction will be weakened compared to the Debye interaction. Judging by Figure 10, the Debye interaction already has a large effect on the diffusion constants for a= 0.5 nm, so we expect further reductions in the ratio of the Lennard-Jones interaction to the Debye interaction to have a limited effect on particle diffusion. Smaller particles should therefore display moderately reduced diffusion constants. A moderate decrease in diffusion constant would lead to a moderately heightened resistance, as seen by Equations (24) and (26).

The size of the Debye interaction depends on the charges involved. CS in the brain tends to be singly sulfated and consequently has twice the charge of HA (Logsdon et al., 2022). One step in moving toward a more detailed model of PNNs could therefore include higher charges on the chain beads, as CSPGs are attached to HA. An increased bead charge will increase the strength of the Debye interaction, as seen by Equation (4). Similarly, the Debye strength was increased in Section 4.7, where the charge of the particle was varied. As seen from Figure 10, the diffusion constant did not depend strongly on particle charge and hence the size of the Debye interaction, though charged particles exhibited stronger diffusion constants than neutral ones. These findings indicate that the exact strength of the Debye interaction is not a strong determinant of diffusion.

However, the polymer chains will experience a stronger interchain repulsion upon an increased bead charge, giving rise to a somewhat more stretched conformation of the chains. As seen in Figure 6, stretching of chains leads to an increase in the diffusion coefficient perpendicular to the wall. For the charged particle, this will decrease the resistance. Figure 6 and Supplementary Figure S6 indicate that the specific geometry of the system plays a relatively small role in diffusion. As the Debye interaction decays quickly with separation, it is fair to assume that the stretching of the chains should be rather moderate. We therefore only expect a slight increase in the diffusion constant and decrease in resistance when the bead charge is increased.

Another simplification made in the model is that only one charged particle diffuses within the brush. The ionic concentration in the brain is on the order of magnitude of 100 mM (Raimondo et al., 2015), meaning that ions present in the fluid interact with each other. This interaction can affect transport, most likely lowering the diffusion constant as the particle is subject to forces from a lot of other particles, effectively slowing it down (Sterratt et al., 2011). The resistance would then be increased. Given the order of magnitude difference between the brush resistance and the membrane resistance of a neuron, we believe that neither of the additional effects discussed should increase the PNN resistance to an extent that it is comparable to the membrane resistance.

6. Conclusion

Diffusion was characterized for a freely diffusing particle in a system of negatively charged polymers tethered to a substrate, a simplified model of the perineuronal nets. The anisotropic nature of the system yielded different diffusion constants in directions parallel and perpendicular to the substrate for denser brushes. Both diffusion constants displayed a clear reduction compared to bulk, particularly for smaller brush spacings. This finding supports the notion that PNNs restrict transport.

Isotropic diffusion was recovered in the limit of very dilute brushes, where the system is close to bulk. In order to characterize how the system approached bulk, the ratio of the diffusion constants were fit to a decaying power law.

To study the effect of brush geometry and dynamics on diffusion, several different system types were probed. Altering the dynamics of the brush yielded only small differences in diffusion of the free particle, with the diffusion constants of static brushes and dynamic brushes with and without a bending term all agreeing within the standard deviation. The free particle diffusion constant in a system of straight, fixed chains yielded a higher diffusion constant in the perpendicular direction than for the other systems, while the parallel diffusion constant remained the same within the standard deviation. In summary, the exact dynamics of the brush did not affect the diffusion constants noticeably.

The behavior of each diffusion constant with increasing spacing d was compared to existing theories for intermediate values of d in a heuristic manner. The hyperbola of revolution and ordered packings models provided the best fit of the theoretical models considered to the diffusion constant as a function of porosity. A modified version of the cation-exchange resin model was implemented and displayed an acceptable performance. A power law proved the best fit to D vs. d in the perpendicular direction for both system types and was used to estimate the diffusion time through a h = 500 nm brush.

Simulations were performed for different particle charges and combined with theory to obtain electrical properties of the brush. The resistance of the brush was found to be orders of magnitude smaller than the resistance of a neuron membrane, implying that the PNNs should not affect the resistance of a cell to a large extent.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10.5281/zenodo.7304745.

Author contributions

KH performed simulations and analysis. AM-S and KH contributed equally to modeling, concepts and writing. All authors contributed to the article and approved the submitted version.

Funding

This study was funded by The Research Council of Norway (https://www.forskningsradet.no/en/) grant no. 568117 to KH and AM-S.

Acknowledgments

We thank Gaute Einevoll and Jonas Verhellen for enlightening discussions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2022.967735/full#supplementary-material

References

Akanni, K. A., Evans, J. W., and Abramson, I. S. (1987). Effective transport coefficients in heterogeneous media. Chem. Eng. Sci. 42, 1945–1954. doi: 10.1016/0009-2509(87)80141-0

CrossRef Full Text | Google Scholar

Attili, S., Borisov, O. V., and Richter, R. P. (2012). Films of end-grafted hyaluronan are a prototype of a brush of a strongly charged, semiflexible polyelectrolyte with intrinsic excluded volume. Biomacromolecules 13, 1466–1477. doi: 10.1021/bm3001759

PubMed Abstract | CrossRef Full Text | Google Scholar

Bano, F., Banerji, S., Howarth, M., Jackson, D. G., and Richter, R. P. (2016). A single molecule assay to probe monovalent and multivalent bonds between hyaluronan and its key leukocyte receptor CD44 under force. Sci. Rep. 6, 34176. doi: 10.1038/srep34176

PubMed Abstract | CrossRef Full Text | Google Scholar

Bathe, M., Rutledge, G. C., Grodzinsky, A. J., and Tidor, B. (2005). A coarse-grained molecular model for glycosaminoglycans: application to chondroitin, chondroitin sulfate, and hyaluronic acid. Biophys. J. 88, 3870–3887. doi: 10.1529/biophysj.104.058800

PubMed Abstract | CrossRef Full Text | Google Scholar

Beeckman, J. W. (1990). Mathematical description of heterogeneous materials. Chem. Eng. Sci. 45, 2603–2610. doi: 10.1016/0009-2509(90)80148-8

CrossRef Full Text | Google Scholar

Cai, L.-H., Panyukov, S., and Rubinstein, M. (2011). Mobility of nonsticky nanoparticles in polymer liquids. Macromolecules 44, 7853–7863. doi: 10.1021/ma201583q

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, Q., Zuo, C., Li, L., and He, H. (2010). Shearing and compression behavior of end-grafted polyelectrolyte brushes with mono- and trivalent counterions: a molecular dynamics simulation. Model. Simulat. Mater. Sci. Eng. 18, 075001. doi: 10.1088/0965-0393/18/7/075001

CrossRef Full Text | Google Scholar

Carbajal-Tinoco, M. D., Lopez-Fernandez, R., and Arauz-Lara, J. L. (2007). Asymmetry in colloidal diffusion near a rigid wall. Phys. Rev. Lett. 99, 138303. doi: 10.1103/PhysRevLett.99.138303

PubMed Abstract | CrossRef Full Text | Google Scholar

Carulli, D., Broersen, R., de Winter, F., Muir, E. M., Meškovi,ć, M., de Waal, M., et al. (2020). Cerebellar plasticity and associative memories are controlled by perineuronal nets. Proc. Natl. Acad. Sci. U.S.A. 117, 6855–6865. doi: 10.1073/pnas.1916163117

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Xu, H., Ma, Y., Liu, J., and Zhang, L. (2022). Diffusion of polymer-grafted nanoparticles with dynamical fluctuations in unentangled polymer melts. Phys. Chem. Chem. Phys. 24, 11322–11335. doi: 10.1039/D2CP00002D

PubMed Abstract | CrossRef Full Text | Google Scholar

Csajka, F. S., and Seidel, C. (2000). Strongly charged polyelectrolyte brushes: a molecular dynamics study. Macromolecules 33, 2728–2739. doi: 10.1021/ma990096l

CrossRef Full Text | Google Scholar

de Beer, S., Mensink, L. I. S., and Kieviet, B. D. (2016). Geometry-dependent insertion forces on particles in swollen polymer brushes. Macromolecules 49, 1070–1078. doi: 10.1021/acs.macromol.5b01960

CrossRef Full Text | Google Scholar

Deepa, S. S., Carulli, D., Galtrey, C., Rhodes, K., Fukuda, J., Mikami, T., et al. (2006). Composition of perineuronal net extracellular matrix in rat brain: a different disaccharide composition for the net-associated proteoglycans. J. Biol. Chem. 281, 17789–17800. doi: 10.1074/jbc.M600544200

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhont, J. K. (1996). An Introduction to Dynamics of Colloids, volume 2 of em Studies in Interface Science. 1st Edn. Amsterdam, Elsevier.

Google Scholar

Doyle, P. S., Shaqfeh, E. S. G., and Gast, A. P. (1998). Rheology of Polymer Brushes: A Brownian 695 Dynamics Study. Macromolecules 31, 5474–5486. doi: 10.1021/ma970821x Publisher: American 696 Chemical Society.

CrossRef Full Text | Google Scholar

Edward, J. T. (1970). Molecular volumes and the stokes-einstein equation. J. Chem. Educ. 47, 261. doi: 10.1021/ed047p261

CrossRef Full Text | Google Scholar

Einstein, A. (1905). On the movement of small particles suspended in stationary liquids required by the molecular-kinetic theory of heat. Ann. Physik 17, 549–560.

Google Scholar

Fawcett, J. W., Oohashi, T., and Pizzorusso, T. (2019). The roles of perineuronal nets and the perinodal extracellular matrix in neuronal function. Nat. Rev. Neurosci. 20, 451–465. doi: 10.1038/s41583-019-0196-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferguson, B. R., and Gao, W.-J. (2018). PV interneurons: critical regulators of E/I balance for prefrontal cortex-dependent behavior and psychiatric disorders. Front. Neural Circ. 12, 37. doi: 10.3389/fncir.2018.00037

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghanbarian, B., Hunt, A. G., Ewing, R. P., and Sahimi, M. (2013). Tortuosity in porous media: a critical review. Soil Sci. Soc. Am. J. 77, 1461–1477. doi: 10.2136/sssaj2012.0435

CrossRef Full Text | Google Scholar

Gheorghiu, S., and Coppens, M.-O. (2004). Heterogeneity explains features of “anomalous” thermodynamics and statistics. Proc. Natl. Acad. Sci. U.S.A. 101, 15852–15856. doi: 10.1073/pnas.0407191101

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., et al. (2020). Array programming with NumPy. Nature 585, 357–362. doi: 10.1038/s41586-020-2649-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hehmeyer, O. J., and Stevens, M. J. (2005). Molecular dynamics simulations of grafted polyelectrolytes on two apposing walls. J. Chem. Phys. 122, 134909. doi: 10.1063/1.1871937

PubMed Abstract | CrossRef Full Text | Google Scholar

Horkay, F., Chremos, A., Douglas, J. F. L., Jones, R., Lou, J., et al. (2020). Systematic investigation of synthetic polyelectrolyte bottlebrush solutions by neutron and dynamic light scattering, osmometry, and molecular dynamics simulation. J. Chem. Phys. 152, 194904. doi: 10.1063/5.0007271

PubMed Abstract | CrossRef Full Text | Google Scholar

Jouhanneau, J.-S., Kremkow, J., and Poulet, J. F. A. (2018). Single synaptic inputs drive high-precision action potentials in parvalbumin expressing GABA-ergic cortical neurons in vivo. Nat. Commun. 9, 1540. doi: 10.1038/s41467-018-03995-2

PubMed Abstract | CrossRef Full Text

Kamerlin, N. (2017). Computer Simulations of Polymer Gels : Structure, Dynamics, and Deformation (Ph.D. thesis). Publisher: Acta Universitatis Upsaliensis.

Google Scholar

Kim, M., Schmitt, S. K., Choi, J. W., Krutty, J. D., and Gopalan, P. (2015). From self-assembled monolayers to coatings: advances in the synthesis and nanobio applications of polymer brushes. Polymers 7, 1346–1378. doi: 10.3390/polym7071346

CrossRef Full Text | Google Scholar

Kinjo, T., Yoshida, H., and Washizu, H. (2018). Coarse-grained simulations of polyelectrolyte brushes using a hybrid model. Colloid Polym. Sci. 296, 441–449. doi: 10.1007/s00396-017-4258-7

CrossRef Full Text | Google Scholar

Kremer, K., and Grest, G. S. (1990). Dynamics of entangled linear polymer melts: a molecular-dynamics simulation. J. Chem. Phys. 92, 5057–5086. doi: 10.1063/1.458541

CrossRef Full Text | Google Scholar

Lemons, D. S., and Gythiel, A. (1997). Paul Langevin's 1908 paper “On the Theory of Brownian Motion” [“Sur la théorie du mouvement brownien,” C. R. Acad. Sci. (Paris) 146, 530–533 (1908)]. Am. J. Phys. 65, 1079–1081. doi: 10.1119/1.18725

CrossRef Full Text | Google Scholar

Li, L., Cao, Q., and Zuo, C. (2018). Effect of counterion valence on conformational behavior of spherical polyelectrolyte brushes confined between two parallel walls. Polymers 10, 363. doi: 10.3390/polym10040363

PubMed Abstract | CrossRef Full Text | Google Scholar

Logsdon, A. F., Francis, K. L., Richardson, N. E., Hu, S. J., Faber, C. L., Phan, B. A., et al. (2022). Decoding perineuronal net glycan sulfation patterns in the Alzheimer's disease brain. Alzheimers Demen. 18, 942–954. doi: 10.1002/alz.12451

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, Y., Wang, C., Pang, A.-P., Zhang, X., Wang, D., and Lu, X. (2021). Low-concentration salt solution changes the interfacial molecular behavior of polyelectrolyte brushes. Macromolecules 54, 6006–6013. doi: 10.1021/acs.macromol.1c00119

CrossRef Full Text | Google Scholar

Mackie, J. S., and Meares, P. (1955). The diffusion of electrolytes in a cation-exchange resin membrane. I. theoretical. Proc. R. Soc. Lond. A Math. Phys. Sci. 232, 498–509. doi: 10.1098/rspa.1955.0234

CrossRef Full Text | Google Scholar

Merlitz, H., Wu, C.-X., and Sommer, J.-U. (2012). Inclusion free energy of nanoparticles in polymer brushes. Macromolecules 45, 8494–8501. doi: 10.1021/ma301781b

PubMed Abstract | CrossRef Full Text | Google Scholar

Morawski, M., Reinert, T., Meyer-Klaucke, W., Wagner, F. E., Tröger, W., Reinert, A., et al. (2015). Ion exchanger in the brain: quantitative analysis of perineuronally fixed anionic binding sites suggests diffusion barriers with ion sorting properties. Sci. Rep. 5, 16471. doi: 10.1038/srep16471

PubMed Abstract | CrossRef Full Text | Google Scholar

Neale, G. H., and Nader, W. K. (1973). Prediction of transport processes within porous media: diffusive flow processes within an homogeneous swarm of spherical particles. AIChE J. 19, 112–119. doi: 10.1002/aic.690190116

CrossRef Full Text | Google Scholar

Nelson, P. (2014). Biological Physics: Energy, Information, Life. New York, NY: W. H. Freedman and Company.

Google Scholar

Ou, Y., Sokoloff, J. B., and Stevens, M. J. (2012). Comparison of the kinetic friction of planar neutral and polyelectrolyte polymer brushes using molecular dynamics simulations. Phys. Rev. E 85, 011801. doi: 10.1103/PhysRevE.85.011801

PubMed Abstract | CrossRef Full Text | Google Scholar

Petersen, E. E. (1958). Diffusion in a pore of varying cross section. AIChE J. 4, 343–345. doi: 10.1002/aic.690040322

CrossRef Full Text | Google Scholar

Pizzorusso, T., Medini, P., Berardi, N., Chierzi, S., Fawcett, J. W., and Maffei, L. (2002). Reactivation of ocular dominance plasticity in the adult visual cortex. Science 298, 1248–1251. doi: 10.1126/science.1072699

PubMed Abstract | CrossRef Full Text | Google Scholar

Postl, A., Hilgert, P. P. P., Markevich, A., Madsen, J., Mustonen, K., Kotakoski, J., et al. (2022). Indirect measurement of the carbon adatom migration barrier on graphene. Carbon 196, 596–601. doi: 10.1016/j.carbon.2022.05.039

CrossRef Full Text | Google Scholar

Qian, N., and Sejnowski, T. J. (1988). “Electrodiffusion model of electrical conduction in neuronal processes,” in Cellular Mechanisms of Conditioning and Behavioral Plasticity, eds C. D. WoodyD. L. Alkon, and J. L. McGaugh (Boston, MA: Springer US), 237–244.

Google Scholar

Quan, T., Zheng, T., Yang, Z., Ding, W., Li, S., Li, J., et al. (2013). NeuroGPS: automated localization of neurons for brain circuits using L1 minimization model. Sci. Rep. 3, 1414. doi: 10.1038/srep01414

PubMed Abstract | CrossRef Full Text | Google Scholar

Raimondo, J. V., Burman, R. J., Katz, A. A., and Akerman, C. J. (2015). Ion dynamics during seizures. Front. Cell. Neurosci. 9, 419. doi: 10.3389/fncel.2015.00419

PubMed Abstract | CrossRef Full Text | Google Scholar

Richter, R. P., Baranova, N. S., Day, A. J., and Kwok, J. C. (2018). Glycosaminoglycans in extracellular matrix organisation: are concepts from soft matter physics key to understanding the formation of perineuronal nets? Curr. Opin. Struct. Biol. 50, 65–74. doi: 10.1016/j.sbi.2017.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, L., and Chen, Z. (2007). Critical review of the impact of tortuosity on diffusion. Chem. Eng. Sci. 62, 3748–3755. doi: 10.1016/j.ces.2007.03.041

CrossRef Full Text | Google Scholar

Singh, M. K., Ilg, P., Espinosa-Marsal, R. M., Kröger, M., and Spencer, N. D. (2015). Polymer brushes under shear: molecular dynamics simulations compared to experiments. Langmuir 31, 4798–4805. doi: 10.1021/acs.langmuir.5b00641

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, M. K., Ilg, P., Espinosa-Marzal, R. M., Spencer, N. D., and Kröger, M. (2016). Influence of chain stiffness, grafting density and normal load on the tribological and structural behavior of polymer brushes: a nonequilibrium-molecular-dynamics study. Polymers 8, 254. doi: 10.3390/polym8070254

PubMed Abstract | CrossRef Full Text | Google Scholar

Spirin, L., Galuschko, A., Kreer, T., Binder, K., and Baschnagel, J. (2011). Polymer-brush lubricated surfaces with colloidal inclusions under shear inversion. Phys. Rev. Lett. 106, 168301. doi: 10.1103/PhysRevLett.106.168301

PubMed Abstract | CrossRef Full Text | Google Scholar

Sterratt, D., Graham, B., Gillies, A., and Willshaw, D. (2011). Principles of Computational Modelling in Neuroscience, 1st Edn. Cambridge: Cambridge University Press.

Google Scholar

Stukowski, A. (2009). Visualization and analysis of atomistic simulation data with OVITO-the open visualization tool. Model. Simulat. Mater. Sci. Eng. 18, 015012. doi: 10.1088/0965-0393/18/1/015012

PubMed Abstract | CrossRef Full Text | Google Scholar

Suttkus, A., Rohn, S., Weigel, S., Glöckner, P., Arendt, T., and Morawski, M. (2014). Aggrecan, link protein and tenascin-R are essential components of the perineuronal net to protect neurons against iron-induced oxidative stress. Cell Death Dis. 5, e1119-e1119. doi: 10.1038/cddis.2014.25

PubMed Abstract | CrossRef Full Text | Google Scholar

Syková, E., and Nicholson, C. (2008). Diffusion in brain extracellular space. Physiol. Rev. 88, 1277–1340. doi: 10.1152/physrev.00027.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Tartakovsky, D. M., and Dentz, M. (2019). Diffusion in porous media: phenomena and mechanisms. Transport Porous Media 130, 105–127. doi: 10.1007/s11242-019-01262-6

CrossRef Full Text | Google Scholar

Thompson, A. P., Aktulga, H. M., Berger, R., Bolintineanu, D. S., Brown, W. M., Crozier, P. S., et al. (2022). LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comput. Phys. Commun. 271, 108171. doi: 10.1016/j.cpc.2021.108171

CrossRef Full Text | Google Scholar

Thompson, E. H., Lensj,ø, K. K., Wigestrand, M. B., Malthe-Sørenssen, A., Hafting, T., and Fyhn, M. (2018). Removal of perineuronal nets disrupts recall of a remote fear memory. Proc. Natl. Acad. Sci. U.S.A. 115, 607–612. doi: 10.1073/pnas.1713530115

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsien, R. Y. (2013). Very long-term memories may be stored in the pattern of holes in the perineuronal net. Proc. Natl. Acad. Sci. U.S.A. 110, 12456–12461. doi: 10.1073/pnas.1310158110

PubMed Abstract | CrossRef Full Text | Google Scholar

van 't Spijker, H. M., and Kwok, J. C. F. (2017). A sweet talk: the molecular systems of perineuronal nets in controlling neuronal communication. Front. Integr. Neurosci. 11, 33. doi: 10.3389/fnint.2017.00033

PubMed Abstract | CrossRef Full Text | Google Scholar

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261–272. doi: 10.1038/s41592-019-0686-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Wetzler, S. P., Miller, K. A., Kisley, L., Stanton, A. L. D., Braun, P. V., and Bailey, R. C. (2020). Real-time measurement of polymer brush dynamics using silicon photonic microring resonators: analyte partitioning and interior brush kinetics. Langmuir 36, 10351–10360. doi: 10.1021/acs.langmuir.0c01336

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Lu, T. J., and Kim, T. (2014). An analytical model for permeability of isotropic porous media. Phys. Lett. A 378, 2308–2311. doi: 10.1016/j.physleta.2014.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Zelenovskiy, P. S., Domingues, E. M., Slabov, V., Kopyl, S., Ugolkov, V. L., Figueiredo, F. M. L., et al. (2020). Efficient water self-diffusion in diphenylalanine peptide nanotubes. ACS Appl. Mater. Interfaces 12, 27485–27492. doi: 10.1021/acsami.0c03658

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Q., and Xiang, X. (2013). Adsorption of a spherical nanoparticle in polymer brushes: brownian dynamics investigation. Physica A 392, 3857–3862. doi: 10.1016/j.physa.2013.05.001

CrossRef Full Text | Google Scholar