# Reaction–Diffusion Model-Based Research on Formation Mechanism of Neuron Dendritic Spine Patterns

^{1}Institute of Robotics and Automatic Information Systems, College of Artificial Intelligence, Nankai University, Tianjin, China^{2}State Key Laboratory of Medicinal Chemical Biology, School of Medicine, Nankai University, Tianjin, China

The pattern abnormalities of dendritic spine, tiny protrusions on neuron dendrites, have been found related to multiple nervous system diseases, such as Parkinson's disease and schizophrenia. The determination of the factors affecting spine patterns is of vital importance to explore the pathogenesis of these diseases, and further, search the treatment method for them. Although the study of dendritic spines is a hot topic in neuroscience in recent years, there is still a lack of systematic study on the formation mechanism of its pattern. This paper provided a reinterpretation of reaction-diffusion model to simulate the formation process of dendritic spine, and further, study the factors affecting spine patterns. First, all four classic shapes of spines, mushroom-type, stubby-type, thin-type, and branched-type were reproduced using the model. We found that the consumption rate of substrates by the cytoskeleton is a key factor to regulate spine shape. Moreover, we found that the density of spines can be regulated by the amount of an exogenous activator and inhibitor, which is in accordance with the anatomical results found in hippocampal CA1 in SD rats with glioma. Further, we analyzed the inner mechanism of the above model parameters regulating the dendritic spine pattern through Turing instability analysis and drew a conclusion that an exogenous inhibitor and activator changes Turing wavelength through which to regulate spine densities. Finally, we discussed the deep regulation mechanisms of several reported regulators of dendritic spine shape and densities based on our simulation results. Our work might evoke attention to the mathematic model-based pathogenesis research for neuron diseases which are related to the dendritic spine pattern abnormalities and spark inspiration in the treatment research for these diseases.

## Introduction

Dendritic spines are tiny protrusions on neuron dendrites which widely exist in the dendrites of higher animals and play an important role in the formation of most excitatory axodendritic synapses (Harris and Kater, 1994). The function of a spine is related to its shape (Kasai et al., 2003; Bourne and Harris, 2007). Traditionally, there are four basic shapes for dendritic spines: thin-type, stubby-type, mushroom-type, and branched-type (González-Tapia et al., 2016; Luczynski et al., 2016). Among them, thin dendritic spines show high plasticity and are related to learning, while mushroom dendritic spines show weak plasticity and are related to memory function. In addition, the density of spines directly influences the density of synapses. Researchers have found that pattern abnormalities of dendritic spine, especially the abnormal proportion of various types of dendritic spines and density variation of dendritic spines, were related to multiple nervous system diseases. For example, Pyronneau et al. reported an overabundance of thin-type spines, a kind of immature dendritic spines, in the somatosensory cortex of Fragile X syndrome model mice (Pyronneau et al., 2017). It has been reported that there are striatal dendrites with few dendritic spines in Parkinson's disease (McNeill et al., 1988). It was also been found that reduced dendritic spine density in individuals with schizophrenia (Glantz and Lewis, 2000; Sweet et al., 2008) and Huntington's disease (Richards et al., 2011). Also, it is recognized that dendritic spine loss is an early feature of Alzheimer's disease (Kommaddi et al., 2018; O'Neal et al., 2018). Thus, the exploration of shape and density factors of dendritic spines is of vital importance to understand the pathogenesis of these diseases, and further, search the treatment method for them.

The current research on dendritic spines pattern is mainly performed by statically observing the cerebral cortex in animals (Kommaddi et al., 2018; Ratliff et al., 2019). It has been confirmed that the pattern of dendritic spines is influenced by neuron activity (Portera-Cailliau et al., 2003; González-Tapia et al., 2016) and some substances, such as drebrin (Hayashi et al., 1996), Rho GTPase Rac1 (Pyronneau et al., 2017) and F-actin (Kommaddi et al., 2018). The above researches usually only proposed one factor of dendritic spine patterns once while the pattern formation of dendritic spines is a dynamic process involving a variety of chemical reactions that are regulated by multiple factors. In summary, there is still a lack of systematic study on the mechanism of pattern formation showing influences of multiple factors on the formed pattern of dendritic spines.

Mathematic modeling on dendritic spines development has become an important tool to study the structure and plasticity of dendritic spines in recent years. For example, Kasai et al. used the volume of dendritic spines as an index to measure the structure of dendritic spines and applied the Brownian motion model to simulate the volume of dendritic spines, exploring the close relationship between spine structure and function (Kasai et al., 2010). The Brownian motion model describes a random phenomenon, but the pattern formation of dendritic spines is a process regulated by gene and environment instead of a random process, making that model unsuitable for simulating the pattern formation. Besides, Miermans et al. simulated dendritic spine membranes during shape alternation using the Canham-Helfrich energy functional, which is used to describe the relationship between the bending rigidity of the membrane and the force generated by the cytoskeleton (Miermans et al., 2017). Their results demonstrate that the cytoskeleton is a key factor in determining the shape of dendritic spines, but this model lacks an explanation for the change in cytoskeletons, and their hypothesis about the approximate rotational symmetry of dendritic spines seems inapplicable to branched-type dendritic spines. Varner et al. explained the process of epithelial cell formation patterns using four mechanisms: cell division, cell insertion, cell deformation, and media filling (Varner and Nelson, 2014). However, these explanations cannot be applied in the study of sub-cellular structures such as dendritic spines.

In Turing theory, if the chemical substances involved in the interaction have diffusion, the original equilibrium state will be broken, which is called Turing instability (Turing, 1952). The reaction-diffusion model (Gierer and Meinhardt, 1972; Meinhardt, 1976), based on Turing's theory, illuminates the reactions between chemical substances in developing biological systems. It has been utilized to simulate Pomacanthus skin stripe patterns (Kondo and Asai, 1995), vascular mesenchymal cells patterns (Garfinkel et al., 2004), mouse limb development (Miura et al., 2006), lung branching patterns (Guo et al., 2014a; Hagiwara et al., 2015), and self-organizing morphogenesis (Okuda et al., 2018; Landge et al., 2020). In our previous work, side branching and tip branching of the lung were investigated using the reaction-diffusion model, which was verified by spatiotemporal parameters (Guo et al., 2014a). However, the patterns developed in previous work were not enough to describe the complex patterns in dendritic spines. Because different from the obtained side branches which were equally spaced, the dendritic spines studied in this paper are usually uneven. In spite of its potential use in simulate branching patterns, the strong non-linearity of the reaction-diffusion model makes it difficult to intuitively draw the relationship between parameter values and simulation results, which is inconvenient for the analysis of the inner mechanism of the model. Addressing this problem, dispersion relation was used to analyze Turing instability (Guo et al., 2014b; Saleem and Ali, 2018) to prove the mathematical mechanism of the simulation results. In previous research, we have investigated the mathematic mechanism through Turing instability analysis and found that different Turing wavelengths are underlying the different patterns in a lung (Xu et al., 2017). However, the relationship between Turing wavelength and branch density has not been investigated yet.

This paper reinterpreted the traditional reaction-diffusion model through the introduction of exogenous activator term and exogenous inhibitor term to simulate the formation process of dendritic spine, and further, study the factors affecting spine patterns. All four spine shapes, mushroom-type, stubby-type, thin-type, and branched-type, were reproduced using the model. Further, we found that the consumption rate of substrates by the cytoskeleton regulates the shape. Secondly, we found that the addition of an exogenous activator causes the spines to become denser, while the addition of an exogenous inhibitor causes the spines to become sparser, which provided a potential explanation for the anatomical results that spine decrease in hippocampal CA1 in SD rats with glioma. Finally, through Turing instability analysis, we found that Turing wavelength variation is the deep mathematical mechanism behind above parameters regulating spine density. Namely, the addition of an exogenous activator decreases the Turing wavelength, causing the density of the dendritic spines to increase, while the addition of an exogenous inhibitor increases the Turing wavelength, causing the density of the dendritic spines to decrease. Finally, the deep regulation mechanisms of several regulators of dendritic spine shape and density reported in other references were discussed based on our simulation results. We hope that our work could evoke attention to the mathematic model-based research for neuron diseases related to the dendritic spine pattern abnormalities and spark inspiration in the treatment research for these diseases.

## Materials and Methods

### Reaction-Diffusion Model

The reaction-diffusion model is defined by Equation (1) (Meinhardt, 1976). It is a group of partial differential equations describing the reactions between activator *A*, inhibitor *H*, substrate *S*, and cytoskeleton *Y*.

The reaction-diffusion model illuminates the reactions between chemical substances in developing biological systems. According to this model, neurons express activators (at a rate ρ_{A}) and inhibitors (at a rate ρ_{H}). Activators behave with self-catalysis (at a rate *c*) and catalyze inhibitors (at a rate *c*), while inhibitors inhibit activators. Simultaneously, activators and inhibitors behave with degradation and diffusion (activators degrade at a rate μ and diffuse at a rate *D*_{A}, whereas inhibitors degrade at a rate υ and diffuse at a rate *D*_{H}). High concentrations of activator accelerate the polymerization of cytoskeletons, inducing the development of dendritic spines. Because the diffusion rate of inhibitors is higher than that of activator, the polymerization of the cytoskeleton in the growth center is accelerated, and the polymerization of the cytoskeleton outside the growing center is inhibited. Thus, the dendritic spine grows in a certain direction, instead of displaying isotropous growth. The neuron creates substrate (at a rate *c*_{0}), while the cytoskeleton consumes substrate (at a rate ε). Substrate accelerates the catalysis of the activator. Similarly, the substrate behaves via degradation and diffusion (degrades at a rate γ and diffuses at a rate *D*_{S}), as well. Because the synthesis of the cytoskeleton consumes substrate, the peak concentration areas of activators and inhibitors, as well as the cytoskeleton, move in the direction of high substrate concentrations (Figure 1).

**Figure 1**. Schematic of the development process. The neuron expresses activators and inhibitors. Activators gather at the tip, while inhibitors diffuse into the surrounding area due to a higher diffusion rate, making only the tip develops. This mechanism makes dendritic spines grow in a certain direction instead of exhibiting isotropous growth.

The development patterns of dendritic spines are determined by the neuron activity (Bloodgood and Sabatini, 2005) and the exogenous substances. The neuron activity is described by the rate of substrate consumed by the cytoskeleton (ε) in our model. Exogenous substances include exogenous activators and exogenous inhibitors (Kommaddi et al., 2018). To describe the influence of exogenous substances, we added the exogenous activator term (δ_{A}) and the exogenous inhibitor term (δ_{H}) into the reaction-diffusion model:

The new model includes 16 parameters, most of which are fixed parameters, such as reaction-term parameter *c*, degradation-term parameters μ, υ, and γ, diffusion-term parameters *D*_{A}, *D*_{H}, and *D*_{S}, and growth-term parameters *d, e*, and *f* . The values of fixed parameters are decided by the chemical characteristics of substances or cells, and the model has been proven to be robust to perturbations of fixed parameters (Murray, 1982). The other parameters are variable (ρ_{A}, δ_{A}, ρ_{H}, δ_{H}, *c*_{0}, and ε), whose values depend on the condition of the development system. In this work, we studied the effect of the neuron activity and the exogenous substances on dendritic spines. Thus, we set parameters δ_{A}, δ_{H}, and ε in Equation (2) as variable parameters.

The values of parameters were set according to previous research in lung branching patterns. In previous work, we set the values of parameters as: *c* = 0.002, μ = 0.18, υ = 0.04, ρ_{A} = 0.063, ρ_{H} = 0.00005, *c*_{0} = 0.02, γ = 0.02, ε = 0.045, *D*_{A} = 0.02, *D*_{H} = 0.32, *D*_{S} = 0.06, *d* = 0.0033, *e* = 0.1, and *f* = 10. We verified the consistency of the mathematical model under certain parameters with the actual biological process by converting the time and space in the numerical simulation and comparing them with the spatiotemporal scale of real lung development (Guo et al., 2014a). The values of fixed parameters and the value ranges of variable parameters in the lung branching model provide references in our new model.

### Numerical Simulation

In this work, we investigated the factors of shape and density of spines using a reaction-diffusion model on different spatial scales. First, we simulated a spine to explore the influence of model parameters on the shape of the spine (Figure 2A). This simulation was performed on a 100 × 100 grid, and the original state was a 10 × 5 pixels rectangular area. Second, we simulated a dendrite with spines to explore the influence of model parameters on the density of spines (Figures 2B,C). This simulation was performed on a 150 × 200 grid, and the original state was a 5 × 10 pixels rectangular area (Figure 2B). Then, a dendrite developed under certain conditions (Figure 2C).

**Figure 2**. The original state of the spine simulation and the dendrite simulation. **(A)** The original state of the spine simulation is used to simulate a single spine in different conditions. Simulations were performed on a 100× 100 grid. The grid size of space is 0.3. Fixed parameters in Equation (2): *c* = 0.002, μ = 0.16, υ = 0.04, ρ_{A} = 0.01, ρ_{H} = 0.00005, *c*_{0} = 0.02, γ = 0.02, *D*_{A} = 0.02, *D*_{H} = 0.26, *D*_{S} = 0.06, *d* = 0.0035, *e* = 0.1, and *f* = 10. **(B)** The first step in the dendrite simulation is used to simulate the dendrite trunk. Simulations were performed on a 150× 200 grid. The grid size of the space is 0.3. Parameters in Equation (2): *c* = 0.002, μ = 0.16, υ = 0.04, ρ_{A} = 0.03, ρ_{H} = 0.0001, δ_{A} = 0, δ_{H} = 0, *c*_{0} = 0.02, ε = 0.017, γ = 0.02, *D*_{A} = 0.02, *D*_{H} = 0.26, *D*_{S} = 0.06, *d* = 0.0035, *e* = 0.1, and *f* = 10. **(C)** The second step in the dendrite simulation grows from **(A)** and is used to simulate spines in different conditions. Fixed parameters in Equation (2): *c* = 0.002, μ = 0.16, υ = 0.04, ρ_{A} = 0.02, ρ_{H} = 0.00005, *c*_{0} = 0.05, γ = 0.02, *D*_{A} = 0.02, *D*_{H} = 0.26, *D*_{S} = 0.06, *d* = 0.0035, *e* = 0.1, and *f* = 10. In **(A,B)**, black regions (*A* = 2, *H* = 0.02, *S* = 1, *Y* = 1) represent a part of a neuron, and white regions (*A* = 0.001, *H* = 0.001, *S* = 1, *Y* = 0) represent the environment surrounding the neuron.

### Turing Instability Analysis Method

To verify the simulation results with mathematics, we explored Turing patterns underlying dendritic spine patterns with our previously developed decoupling method (Guo et al., 2014b). The substrate and cytoskeleton are considered dependent variables of time and space, written as *S*(x, y, t) and *Y*(x, y, t). Then, we put these variables into Equation (2) as parameters and obtained the model of an activator-inhibitor system as:

Branching is a system that can grow and form stable mode, which corresponds to the damped oscillation system of mathematical model. Some points in S-Y space correspond to damped oscillatory systems. The set of these points is called Turing instability space, and the wavelength of damped oscillation system is called Turing wavelength (Turing, 1952). According to its definition, the mathematical expression of Turing space can be calculated. The detailed derivation process is in our previous work (Xu et al., 2017).

To explore dendritic spine development patterns according to Turing instability, a scheme was performed according to the following steps.

• Choose an interesting point (the branching point in the branched spine or a random point on the central axis in others) in a simulation result and plot the *S*-*Y* curve of this point.

• Calculate the Turing instability space using Equation (3).

• Find the intersection of the *S*-*Y* curve and Turing instability space.

• According to the form of the solution of Equation (2), we have

and calculate the dispersion relation:

• Record the maximum of the real part of the eigenvalue (λ_{m}) and corresponding wavenumbers (*k*_{m}).

• Calculate Turing wavelength (Λ) of the point in Step 1:

We used Turing instability analysis to explore the difference of mathematical mechanism behind different patterns of dendritic spines in section Turing Instability Underlying Dendritic Spines.

### Anatomy of Hippocampal CA1 in SD Rat

In this study, images from Golgi-Cox-stained brain slices from SD rats were compared with simulation results. Golgi-Cox staining was carried out with a commercial Golgi staining kit (Keyijiaxin, Tianjin, China). As soon as they were taken from the skulls, the brains were stored in Golgi-Cox staining solution in a dark place for 2 weeks, and the solution was replaced at intervals of 48 h. Then, brain slices were produced using a vibratome (VT 1000S, Leica, Germany) with a thickness of 150 μm. The slices were placed on slides covered with 2% gelatine. Next, the slices were dyed with ammonia for 60 min; washed with water three times; fixed with Kodak film for 30 min; and then washed with water, dehydrated, cleared, and mounted. Later, dendritic spines in the CA1 region of the hippocampus were imaged under the 100 × objective lens with a digital camera. Dendritic trees were detected along CA1 tertiary dendrites derived from secondary dendrites, which started at their point on the primary dendrite. All animal experiments were approved by the Animal Research Ethics Committee, School of Medicine, Nankai University and were performed in accordance with the Animal Management Rules of the Ministry of Health of the People's Republic of China.

## Results

### Dendritic Spine Shape Factors Research Based on Reaction-Diffusion Model

There are four traditional types of dendritic spines: mushroom-type, stubby-type, thin-type, and branched-type (González-Tapia et al., 2016; Luczynski et al., 2016). In order to research the factors of dendritic spine shape, we firstly proposed a classification method of spine shape based on real spine microimages. Then, we classified a spine simulated by our reaction-diffusion model and found the change rule of dendritic spine shape in different conditions.

#### Classification Method of Dendritic Spine Shape

At present, the classification methods of dendritic spines shape are qualitative, expert experience-required. To study the shape of dendritic spines quantitatively, metrics to classify dendritic spines need to be determined. Given a branched-type dendritic spine is easy to identify, here we only propose a classification method for the three types of non-branched spines. First, we measured three geometric qualities of a dendritic spine, namely, the height (*h*), the extreme width of the head (*w*_{head}), and the extreme width of the neck (*w*_{neck}), as shown in Figures 3A–D. Then, based on these three values, we constructed two following dimensionless metrics:

• Relative average width (RAW) measures the overall thickness of spines, defined as

• Relative constriction width (RCW) measures the difference between the head width and the neck width, defined as

We calculated the RAWs and RCWs of eight dendritic spines (including three mushroom-type spines, three stubby-type spines, and two thin-type spines, shown in Figure 3E). Thin-type spines have a thin head and neck, so the value of RAW is small. Both the head and neck of stubby-type spines are thick, and the head is thinner or slightly thicker than the neck, so for them, the value of RAW is usually large, and the value of RCW is small or even negative. For the mushroom-type spines, usually have a large head and a thin neck, their values of RAW and RCW are both large. Based on the above analysis, we set the metrics for three types of dendritic spines. As shown in Figure 3F, the shape differences among the three types of spines are obvious. We chose RAW = 0.4 and RCW = 0.25 as two criteria to classify the three types.

**Figure 3**. Metrics of dendritic spine shape. **(A–D)** Three geometric qualities of dendritic spines, namely, the height (*h*), the extreme width of the head (*w*_{head}), and the extreme width of the neck (*w*_{neck}). **(A)** is a mushroom-type spine, **(B)** is a stubby-type spine, **(C)** is a thin-type spine, and **(D)** is a branched-type spine. For a branched-type spine, the extreme width of the head is meaningless. **(E)** Top: Golgi-Cox staining of brain slices from SD rat hippocampal CA1. Four types of spines emerge in the images. Bottom: We found nine dendritic spines (including three mushroom-type spines, three stubby-type spines, two thin-type spines, and one branched-type spine) in the above images. **(F)** Two-metrics distribution of mushroom-type, stubby-type, and thin-type spines. Three types can be classified by the criteria RAW = 0.4 and RCW = 0.25. **(G)** Flow chart of the classification method of dendritic spine shape.

Finally, we presented a flow chart to distinguish the shapes of dendritic spines (Figure 3G). First, if the dendritic spine has a branching structure, it is recognized as a branched-type spine. Second, if the RAW value is lower than 0.4, it is regarded as a thin-type spine. Finally, if the RCW value is lower than 0.25, it is identified as a stubby-type spine. Otherwise, it is recognized as a mushroom-type spine.

#### Consumption Rate of Substrate Dominates the Spine Shape Based on the Reaction-Diffusion Model

In our previous simulation, the rate that substrates are consumed by cells has been shown to play an important role in the branching pattern (Xu et al., 2017). Thus, we assumed that the consumption rate of substrates, namely, the neuron activity, has an effect on the spine shape. To verify this assumption, we performed the following single-spine simulations.

First, to investigate the influence of the consumption rate of substrates (ε) on the shape of dendritic spine, we adjusted the value of parameter ε in Equation (2). We varied the value of ε from 0.01 to 0.9, and part of the obtained results are shown in Figure 4A (with different amplification factors) (also see Supplementary Videos 1–4, respectively). As the value of ε increases, both RAW and RCW values decrease, and the dendritic spine shapes sequentially undergoes mushroom (0 < ε ≤ 0.02), stubby (0.02 < ε ≤ 0.04), thin (0.04 < ε ≤ 0.7), and branched (0.7 < ε ≤ 0.9) forms. All four dendritic spine shapes can be obtained with an increase in the consumption rate of substrates. This result indicated that neuron activity regulates the shape of dendritic spine.

**Figure 4**. Influence of parameters on the shape of dendritic spines. **(A)** Some of the results of spine simulations for the condition of δ_{A} = 0.01, δ_{H} = 0.00005, 0.01 ≤ ε ≤ 0.9 (with different amplification factors). Arrows mark the newborn parts of spines growing from the original state in Figure 2C. With the enhancement of neuron activity, the values of RAW and RCW decrease, and all four shapes are obtained. The shape can be mushroom (0 < ε ≤ 0.02), stubby (0.02 < ε ≤ 0.04), thin (0.04 < ε ≤ 0.7) or branched (0.7 < ε ≤ 0.9) types. **(B)** The results of spine simulation in the condition of δ_{A} = 0, 0.005, 0.01, 0.015, and 0.02, δ_{H} = 0.00005, ε = 0.03, 0.05, 0.7, and 0.9. With the addition of an exogenous activator, the stubby-type spine becomes mushroom-type, and the thin-type spine becomes stubby-type. **(C)** The results of spine simulation in the condition of δ_{A} = 0.01, δ_{H} = 0, 0.000025, 0.00005, 0.000075, and 0.0001, ε = 0.01, 0.03, 0.07, and 0.7. With the addition of an exogenous inhibitors, the branched-type spine becomes thin-type. In **(A–C)**, fixed parameters in Equation (2): *c* = 0.002, μ = 0.16, υ = 0.04, ρ_{A} = 0.01, ρ_{H} = 0.00005, *c*_{0} = 0.02, γ = 0.02, *D*_{A} = 0.02, *D*_{H} = 0.26, *D*_{S} = 0.06, *d* = 0.0035, *e* = 0.1, and *f* = 10.

In addition, to investigate the influence of exogenous activator (δ_{A}) and exogenous inhibitor (δ_{H}) on the shape of dendritic spine, we adjusted the value of parameter δ_{A} and δ_{H} in Equation (2), respectively. We varied the values of δ_{A} under the conditions of ε = 0.03, 0.05, 0.7, and 0.9 and the values of δ_{H} under the conditions of ε = 0.01, 0.03, 0.07, and 0.7, and the results are shown in Figures 4B,C (also see Supplementary Figures 3, 4, respectively). According to the results, we found that a stubby-type spine transforms to mushroom-type and a thin-type spine transforms to stubby-type with an increase in δ_{A}; additionally, a branched-type spine becomes thin-type with an increase in δ_{H}. However, there is no effect of δ_{A} on branched-type spines and no effect of δ_{H} on mushroom-type and stubby-type spines. These results indicated that both δ_{A} and δ_{H} also regulate the spine shape but they are not dominating factors compared to the consumption rate of substrates.

Therefore, dendritic spines sequentially undergo in-turn transformation of mushroom-type, stubby-type, thin-type, and branched-type, with an increase in the consumption rate of substrates. In contrast, exogenous activators affect non-branched dendritic spines, and exogenous inhibitors affect branched dendritic spines. Thus, the consumption rate of substrates (neuron activity) determines the shape of dendritic spines.

### Dendritic Spine Densities Factors Research Based on Reaction-Diffusion Model

Dendritic spines participate in the formation of most excitatory axodendritic synapses, so the density of spines directly influences the density of synapses. In order to research the factors of dendritic spine density, we simulated a dendrite with spines using the reaction-diffusion model and found the relationship between dendritic spine density and key factors. Moreover, we observed the decrease of spine density in the hippocampal CA1 in rats with glioma and proposed a potential reason for this phenomenon by comparing the simulation results and observation results. Further, we used Turing instability to explain the mathematical mechanism behind the above parameters regulating spine density and found that an exogenous inhibitor and activator changes Turing wavelength through which to regulate spine densities.

#### Exogenous Substances Regulate Spine Density

To investigate the factors of dendritic spine density, we next simulated different spine densities which seen across multiple spines through dendrite simulations. In our previous research, we found that the rates of activator and inhibitor secretion from cells have been shown to play an important role in the density of side branching (Guo et al., 2014a). Similarly, it is reasonable for us to assume that exogenous activator and inhibitor are two key factors influencing the density of dendritic spines.

Firstly, in order to find out the effect of exogenous activator and inhibitor on the spine density, we adjusted the values of the two parameters δ_{A} and δ_{H} based on standard values of δ_{A} = 0.01, δ_{H} = 0.00005, and ε = 1, and we obtained two groups of results (Figures 5A,B). The results showed the density of dendritic spines is positively correlated with δ_{A} and negatively correlated with δ_{H}.

**Figure 5**. Influence of parameters on the density of dendritic spines. **(A)** The results of dendrite simulation under the conditions of δ_{A} = 0, 0.01, 0.02, 0.03, and 0.04, δ_{H} = 0.00005, ε = 1. With the addition of an exogenous activator, the number of dendritic spines increases dramatically. **(B)** The results of dendrite simulation under the conditions of δ_{A} = 0.01, δ_{H} = 0, 0.00005, 0.0001, 0.00015, and 0.0002, ε = 1. With the addition of an exogenous inhibitor, the number of dendritic spines decreases dramatically. **(C)** The results of dendrite simulation under the conditions of δ_{A} = 0.01, δ_{H} = 0.00005, ε = 0.5, 1, 1.5, 2, and 2.5. With increased neuron activity, the shape varies from non-branched-dominant to branched-dominant. Branched-type spines take up more space; thus, the growth of surrounding spines is inhibited. In addition, the parameter ε has no effect on the density; for example, the densities are the same under the conditions of ε = 2 and ε = 2.5. In **(A–C)**, fixed parameters in Equation (2): *c* = 0.002, μ = 0.16, υ = 0.04, ρ_{A} = 0.02, ρ_{H} = 0.00005, *c*_{0} = 0.05, γ = 0.02, *D*_{A} = 0.02, *D*_{H} = 0.26, *D*_{S} = 0.06, *d* = 0.0035, *e* = 0.1, and *f* = 10.

Next, we adjusted the values of the parameter ε to find whether the consumption rate of substrates is another factor of density, and the results are shown in Figure 5C (also see Supplementary Figures 1, 2, respectively). We noticed that the spine shape varied from non-branched-dominant to branched-dominant when ε varies from 0.5 to 2.0. Meanwhile, the spine densities have not significant changes when ε varies; for example, the densities are the same under the conditions of ε = 2 and ε = 2.5.

Through dendrite simulations, we found that exogenous activators increase the density of spines, while exogenous inhibitors have the opposite effect. In comparison to exogenous substances, neuron activity has no significant effect on the density.

#### Application in the Hippocampal CA1 of Rats

The hippocampus plays an important role in memory function and cognitive abilities (Muller et al., 1996). Certain diseases, such as glioma, affect the developmental patterns of dendritic spines on hippocampal neurons. It has also been reported that the impairment of neurocognitive function is a common consequence of glioma, in both glioma patients (Wefel et al., 2016; Van Kessel et al., 2017) and glioma animal models (Wang et al., 2010; Hao et al., 2018). Through anatomy and neuron microimaging (see section Anatomy of Hippocampal CA1 in SD Rat for detail), we found that dendritic spines in rats with glioma were less dense (Figures 6A,B, also see Supplementary Videos 5, 6, respectively).

**Figure 6**. Exogenous inhibitor decreases the spine density in hippocampal CA1 in SD rats with glioma. **(A)** A microscopic image of dendritic spines in the brain of a rat in the sham group. There are dense spines in this image, similar to those in our dendrite simulation results under the condition of δ_{A} = 0.01, δ_{H} = 0, ε = 1. **(B)** A microscopic image of dendritic spines in the brain of a rat in the glioma group. There are sparse spines in this image, similar to those in our dendrite simulation results under the condition of δ_{A} = 0.01, δ_{H} = 0.0001, ε = 1. In **(A,B)**, fixed parameters in Equation (2): *c* = 0.002, μ = 0.16, υ = 0.04, ρ_{A} = 0.02, ρ_{H} = 0.00005, *c*_{0} = 0.05, γ = 0.02, *D*_{A} = 0.02, *D*_{H} = 0.26, *D*_{S} = 0.06, *d* = 0.0035, *e* = 0.1, and *f* = 10.

To study the reasons for various densities in the rat hippocampal CA1, we compared the microscopic images of neurons with our simulation results. It seems that the spine patterns in the brains of the rat sham group were similar to those in the simulation results under the condition of δ_{A} = 0.01, δ_{H} = 0, ε = 1, while the spine patterns in the brains of the rat glioma group were similar to those in the simulation results under the condition of δ_{A} = 0.01, δ_{H} = 0.0001, ε = 1 (Figure 5B). Thus, we considered that the addition of exogenous inhibitors is a potential reason for the decrease of dendritic spine density caused by glioma.

#### Turing Instability Underlying Dendritic Spines

Turing pointed out that the diffusion of chemical substances will break the original equilibrium state of substance concentration, which is called Turing instability (Turing, 1952). Branching patterns can only be generated from models with Turing instability. In order to qualitatively analyze the Turing instability, equilibrium position, and periodicity of the model solution, we have proposed a Turing instability analysis method using dispersion relation in previous research (Guo et al., 2014b) and found that the Turing wavelength is the internal factor causing the change of branching pattern of a lung (Xu et al., 2017) (see section Turing Instability Analysis Method for more details).

Exogenous substances have an effect on the Turing instability space in which a stable pattern can appear, and have no effect on the *S*-*Y* curve that shows the concentration relationship between substrate and cytoskeleton during development, which can be derived from Equation (3). We adjusted the values of the parameters δ_{A} from 0 to 0.2 and then drew an *S*-*Y* curve and Turing instability space in the *S*-*Y* space (Figure 7A). Three intersection points of Turing instability spaces and the corresponding *S*-*Y* curves were marked with black points. These points were substituted into the Equations (4–6) in order to calculate the Turing wavelengths (Figure 7B). An increase in parameter δ_{A} decreases the Turing wavelength. Similarly, we found the intersection points and calculated Turing wavelength under the conditions of δ_{H} = 0, 0.00005, and 0.0001 (Figures 7C,D). An increase in parameter δ_{H} increases the Turing wavelength.

**Figure 7**. Exogenous substances regulate dendritic spine density through Turing wavelength. **(A)** The S-Y curve and Turing instability space in three conditions of parameter δ_{A} (δ_{H} = 0.00005, ε = 1). The yellow, orange, and red regions represent Turing instability spaces in the conditions of δ_{A} = 0, 0.01, and 0.02, respectively. The cyan lines represent S-Y curves, and the black points represent the intersections. An exogenous activator affects the size and location of the Turing instability space, rather than the S-Y curves (artificial errors in choosing intersections result in non-coincidence). **(B)** The dispersion relationship curves in the conditions of δ_{A} = 0, 0.01, and 0.02 (δ_{H} = 0.00005, ε = 1). Adding an exogenous activator decreases the Turing wavelength. **(C)** The S-Y curve and Turing instability space in three conditions of parameter δ_{H} (δ_{A} = 0.01, ε = 1). The yellow, orange, and red regions represent Turing instability spaces in the conditions of δ_{H} = 0, 0.00005, and 0.0001, respectively. The cyan lines represent S-Y curves, and the black points represent the intersections. An exogenous inhibitor affects the size and location of the Turing instability space, rather than the S-Y curves (artificial errors in choosing intersections result in non-coincidence). **(D)** The dispersion relationship curves in the conditions of δ_{H} = 0, 0.00005, and 0.0001 (δ_{A} = 0.01, ε = 1). Adding an exogenous inhibitor increases the Turing wavelength.

As the Turing wavelength implies the spatial periodicity of spines, it is negatively correlated to the density of dendritic spines. In conclusion, exogenous activators make the Turing wavelength smaller and cause an increase in density of dendritic spines, while exogenous inhibitors increase the Turing wavelength and cause a decrease in the density of dendritic spines.

## Discussion

In recent years, various chemicals have been reported to be capable of regulating the process of dendritic spine development. Our research may explore their regulation mechanism in a mathematical view. For example, the actin filaments (F-actin) were considered to be key in regulating the shape of dendritic spines (Miermans et al., 2017). We found the cytoskeleton was one key factor to regulate cell morphology. Hence, F-actin might be considered as the cytoskeleton (*Y*) in our model. It has been found that drebrin is an actin-binding protein in the dendritic spine, and its overexpression causes spine elongation (Hayashi and Shirao, 1999; Koganezawa et al., 2017; Hanamura et al., 2018). Bernstein reported that cofilin severs F-actin, contributing to actin dynamics (Bernstein and Bamburg, 2010). In addition, Calabrese suggested that dendritic spine growth correlates with decreased cofilin activity (Calabrese et al., 2014). According to our simulation results, drebrin and cofilin are similar to the functions of the activator (*A*) and the inhibitor (*H*) in our model respectively. Adenosine-triphosphate (ATP) is closely related to F-actin polymerization and depolymerization (Katkar et al., 2018; Merino et al., 2018), which implies that ATP may correspond to the substrate (*S*) in our model. Based on these hypotheses, we described our inferences as follows: (1) the overexpression of drebrin promotes the binding of F-actin and increases the density of dendritic spines, (2) the overexpression of cofilin hinders the binding of F-actin and decreases the density of dendritic spines, (3) the increase in ATP consumption during the process of creating F-actin results in a different F-actin pattern and causes spines to become mushroom-type, stubby-type, thin-type and branched-type, in turn.

The verification experiments of morphogens is helpful to the correction of model parameters and the support of the conclusion in this work. Here, we proposed two ideas to verify the morphogens mentioned above: (1) research on the quantitative relationship between spine density and the addition of a substance that influences the expression of drebrin or cofilin, and (2) research on the quantitative relationship between spine shape distribution and ATP consumption during the process of creating F-actin. Moreover, in order to compare the spatiotemporal parameters between simulations and verification experiments quantitatively, 3D simulation is necessary.

With our method, certain diseases could be systematically investigated at the level of chemical reactions. For example, the anomalous rise of rho GTPase Rac1 activity inhibited cofilin in mice with Fragile X syndrome because of a trinucleotide expansion in the FMR1 gene on the X chromosome (Pyronneau et al., 2017). In our model, the decrease in δ_{H} decreases the concentration of the inhibitor (*H*), which results in dense dendritic spines. In another study, the intrathecal administration of latrunculin A, an actin-depolymerizing agent, in mice resulted in a decrease in F-actin levels and symptoms of Alzheimer's disease. Conversely, the intrathecal administration of jasplakinolide, a molecule that stabilizes F-actin, in mice restored F-actin levels and improved symptoms (Kommaddi et al., 2018). The effects of latrunculin A and jasplakinolide are similar to those of exogenous inhibitors and exogenous activators in this model, respectively. Exogenous activators promote the synthesis of the cytoskeleton, while exogenous inhibitors promote the decomposition of the cytoskeleton.

In conclusion, we were devoted to revealing the mechanism of the development patterns of dendritic spines. The results show that the consumption rate of substrate dominates the shape, while the addition of exogenous activators and exogenous inhibitors dominates the density. Our work provided a potential explanation for the phenomenon that sparser spines in the brains of SD rats with glioma and maybe also explain some diseases reported in the literature, such as Fragile X syndrome and Alzheimer's disease. Our research provides novel and fresh insight into the development patterns of dendritic spines, helping search treatment methods for related diseases.

## Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

## Ethics Statement

The animal study was reviewed and approved by Animal Research Ethics Committee, School of Medicine, Nankai University.

## Author Contributions

YJ: conceptualization, methodology, software, formal analysis, investigation, writing-original draft, writing-review and editing, and visualization. QZ: conceptualization, methodology, writing-original draft, writing-review and editing, and funding acquisition. HY: methodology, validation, formal analysis, investigation, writing-original draft, writing-review and editing, and visualization. SG: methodology, software, formal analysis, data curation, and writing-review and editing. MS: validation, formal analysis, writing-review and editing, and funding acquisition. ZY: conceptualization, writing-review and editing, and supervision. XZ: conceptualization, writing-review and editing, supervision, project administration, and funding acquisition. All authors contributed to the article and approved the submitted version.

## Funding

This work was supported by National Natural Science Foundation of China (62027812, U1813210, 61903201) and Opening Foundation of Beijing Advanced Innovation Center for Intelligent Robots and Systems (2019IRS05).

## 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.

## Supplementary Material

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

## References

Bernstein, B. W., and Bamburg, J. R. (2010). ADF/Cofilin: a functional node in cell biology. *Trends Cell Biol.* 20, 187–195. doi: 10.1016/j.tcb.2010.01.001

Bloodgood, B., and Sabatini, B. (2005). Neuronal activity regulates diffusion across the neck of dendritic spines. *Science* 310, 866–869. doi: 10.1126/science.1114816

Bourne, J., and Harris, K. (2007). Do thin spines learn to be mushroom spines? *Curr. Opin. Neurobiol.* 17, 381–386. doi: 10.1016/j.conb.2007.04.009

Calabrese, B., Saffin, J.-M., and Halpain, S. (2014). Activity-dependent dendritic spine shrinkage and growth involve downregulation of cofilin via distinct mechanisms. *PLoS ONE* 9:e94787. doi: 10.1371/journal.pone.0094787

Garfinkel, A., Tintut, Y., Petrasek, D., Boström, K., and Demer, L. L. (2004). Pattern formation by vascular mesenchymal cells. *Proc. Natl. Acad. Sci. U.S.A.* 101, 9247–9250. doi: 10.1073/pnas.0308436101

Gierer, A., and Meinhardt, H. (1972). A theory of biological pattern formation. *Kybernetik* 12, 30–39. doi: 10.1007/BF00289234

Glantz, L. A., and Lewis, D. A. (2000). Decreased dendritic spine density on prefrontal cortical pyramidal neurons in schizophrenia. *JAMA Psychiatry* 57, 65–73. doi: 10.1001/archpsyc.57.1.65

González-Tapia, D., Martínez-Torres, N. I., Hernández-González, M., Guevara, M. A., and González-Burgos, I. (2016). Plastic changes to dendritic spines on layer V pyramidal neurons are involved in the rectifying role of the prefrontal cortex during the fast period of motor learning. *Behav. Brain Res.* 298, 261–267. doi: 10.1016/j.bbr.2015.11.013

Guo, Y., Chen, T.-H., Zeng, X., Warburton, D., Boström, K. I., Ho, C.-M., et al. (2014a). Branching patterns emerge in a mathematical model of the dynamics of lung development. *J. Physiol.* 592, 313–324. doi: 10.1113/jphysiol.2013.261099

Guo, Y., Sun, M., Garfinkel, A., and Zhao, X. (2014b). Mechanisms of side branching and tip splitting in a model of branching morphogenesis. *PLoS ONE* 9:e102718. doi: 10.1371/journal.pone.0102718

Hagiwara, M., Peng, F., and Ho, C.-M. (2015). *In vitro* reconstruction of branched tubular structures from lung epithelial cells in high cell concentration gradient environment. *Sci. Rep.* 5:8054. doi: 10.1038/srep08054

Hanamura, K., Kamata, Y., Yamazaki, H., Kojima, N., and Shirao, T. (2018). Isoform-dependent regulation of drebrin dynamics in dendritic spines. *Neuroscience* 379, 67–76. doi: 10.1016/j.neuroscience.2018.02.038

Hao, S., Gao, J., Hui, W., Zhang, Y., Pavlov, A., Ge, H., et al. (2018). AG-1031 and AG-1503 improve cognitive deficits by promoting apoptosis and inhibiting autophagy in C6 glioma model rats. *Brain Res.* 1699, 1–8. doi: 10.1016/j.brainres.2018.06.026

Harris, K. M., and Kater, S. B. (1994). Dendritic spines: cellular specializations imparting both stability and flexibility to synaptic function. *Annu. Rev. Neurosci.* 17, 341–371. doi: 10.1146/annurev.ne.17.030194.002013

Hayashi, K., Ishikawa, R., Ye, L.-H., He, X.-L., Takata, K., Kohama, K., et al. (1996). Modulatory role of drebrin on the cytoskeleton within dendritic spines in the rat cerebral cortex. *J. Neurosci.* 16, 7161–7170. doi: 10.1523/JNEUROSCI.16-22-07161.1996

Hayashi, K., and Shirao, T. (1999). Change in the shape of dendritic spines caused by overexpression of drebrin in cultured cortical neurons. *J. Neurosci.* 19, 3918–3925. doi: 10.1523/JNEUROSCI.19-10-03918.1999

Kasai, H., Fukuda, M., Watanabe, S., Hayashi-Takagi, A., and Noguchi, J. (2010). Structural dynamics of dendritic spines in memory and cognition. *Trends Neurosci.* 33, 121–129. doi: 10.1016/j.tins.2010.01.001

Kasai, H., Matsuzaki, M., Noguchi, J., Yasumatsu, N., and Nakahara, H. (2003). Structure–stability–function relationships of dendritic spines. *Trends Neurosci.* 26, 360–368. doi: 10.1016/S0166-2236(03)00162-0

Katkar, H. H., Davtyan, A., Durumeric, A. E., Hocky, G. M., Schramm, A. C., Cruz, E. M. D. L., et al. (2018). Insights into the cooperative nature of ATP hydrolysis in actin filaments. *Biophys. J.* 115, 1589–1602. doi: 10.1016/j.bpj.2018.08.034

Koganezawa, N., Hanamura, K., Sekino, Y., and Shirao, T. (2017). The role of drebrin in dendritic spines. *Mol. Cell. Neurosci.* 84, 85–92. doi: 10.1016/j.mcn.2017.01.004

Kommaddi, R. P., Das, D., Karunakaran, S., Nanguneri, S., Bapat, D., Ray, A., et al. (2018). Aβ mediates F-actin disassembly in dendritic spines leading to cognitive deficits in Alzheimer's disease. *J. Neurosci.* 38, 1085–1099. doi: 10.1523/JNEUROSCI.2127-17.2017

Kondo, S., and Asai, R. (1995). A reaction–diffusion wave on the skin of the marine angelfish Pomacanthus. *Nature* 376, 765–768. doi: 10.1038/376765a0

Landge, A., Jordan, B., Diego, X., and Mueller, P. (2020). Pattern formation mechanisms of self-organizing reaction-diffusion systems. *Dev. Biol.* 460, 2–11. doi: 10.1016/j.ydbio.2019.10.031

Luczynski, P., Whelan, S. O., O'Sullivan, C., Clarke, G., Shanahan, F., Dinan, T. G., et al. (2016). Adult microbiota-deficient mice have distinct dendritic morphological changes: differential effects in the amygdala and hippocampus. *Eur. J. Neurosci.* 44, 2654–2666. doi: 10.1111/ejn.13291

McNeill, T. H., Brown, S. A., Rafols, J. A., and Shoulson, I. (1988). Atrophy of medium spiny I striatal dendrites in advanced Parkinson's disease. *Brain Res.* 455, 148–152. doi: 10.1016/0006-8993(88)90124-2

Meinhardt, H. (1976). Models and hypotheses. *Differentiation* 6, 117–123. doi: 10.1111/j.1432-0436.1976.tb01478.x

Merino, F., Pospich, S., Funk, J., Wagner, T., Küllmer, F., Arndt, H.-D., et al. (2018). Structural transitions of F-actin upon ATP hydrolysis at near-atomic resolution revealed by cryo-EM. *Nat. Struct. Mol. Biol.* 25, 528–537. doi: 10.1038/s41594-018-0074-0

Miermans, C. A., Kusters, R. P. T., Hoogenraad, C. C., and Storm, C. (2017). Biophysical model of therole of actin remodeling on dendritic spine morphology. *PLoS ONE* 12:e0170113. doi: 10.1371/journal.pone.0170113

Miura, T., Shiota, K., Morriss-Kay, G., and Maini, P. (2006). Mixed-mode pattern in doublefoot mutant mouse limb–turing reaction-diffusion model on a growing domain during limb development. *J. Theor. Biol.* 240, 562–573. doi: 10.1016/j.jtbi.2005.10.016

Muller, R., Stead, M., and Pach, J. (1996). The hippocampus as a cognitive graph. *J. Gen. Physiol.* 107, 663–694 doi: 10.1085/jgp.107.6.663

Murray, J. (1982). Parameter space for turing instability in reaction diffusion mechanisms: a comparison of models. *J. Theor. Biol.* 98, 143–163. doi: 10.1016/0022-5193(82)90063-7

Okuda, S., Miura, T., Inoue, Y., Adachi, T., and Eiraku, M. (2018). Combining turing and 3d vertex models reproduces autonomous multicellular morphogenesis with undulation, tubulation, and branching. *Sci. Rep.* 8:2386. doi: 10.1038/s41598-018-24858-2

O'Neal, M., Stallings, N., and Malter, J. (2018). Alzheimer's disease, dendritic spines, and calcineurin inhibitors: a new approach? *ACS Chem. Neurosci.* 9, 1233–1234. doi: 10.1021/acschemneuro.8b00213

Portera-Cailliau, C., Pan, D., and Yuste, R. (2003). Activity-regulated dynamic behavior of early dendritic protrusions: evidence for different types of dendritic filopodia. *J. Neurosci. Off. J. Soc. Neurosci.* 23, 7129–7142. doi: 10.1523/JNEUROSCI.23-18-07129.2003

Pyronneau, A., He, Q., Hwang, J.-Y., Porch, M., Contractor, A., and Zukin, R. S. (2017). Aberrant Rac1-cofilin signaling mediates defects in dendritic spines, synaptic function, and sensory perception in fragile X syndrome. *Science Signal.* 10:eaan0852. doi: 10.1126/scisignal.aan0852

Ratliff, W. A., Mervis, R. F., Citron, B. A., Schwartz, B., Rubovitch, V., Schreiber, S., et al. (2019). Mild blast-related TBI in a mouse model alters amygdalar neurostructure and circuitry. *Exp. Neurol.* 315, 9–14. doi: 10.1016/j.expneurol.2019.01.020

Richards, P., Didszun, C., Campesan, S., Simpson, A., Horley, B., Young, K., et al. (2011). Dendritic spine loss and neurodegeneration is rescued by rab11 in models of huntington's disease. *Cell Death Differ.* 18, 191–200. doi: 10.1038/cdd.2010.127

Saleem, M., and Ali, I. (2018). Numerical simulations of turing patterns in a reaction-diffusion model with the chebyshev spectral method. *Eur. Phys. J. Plus* 133:399. doi: 10.1140/epjp/i2018-12265-9

Sweet, R. A., Henteleff, R. A., Zhang, W., Sampson, A. R., and Lewis, D. A. (2008). Reduced dendritic spine density in auditory cortex of subjects with schizophrenia. *Neuropsychopharmacology* 34, 374–389. doi: 10.1038/npp.2008.67

Turing, A. M. (1952). The chemical basis of morphogenesis. *Philos. Trans. R. Soc. Lond.* 237, 37–72. doi: 10.1098/rstb.1952.0012

Van Kessel, E. E., Baumfalk, A. J. E., van Zandvoort, M., Robe, P., and Snijders, T. J. (2017). Tumor-related neurocognitive dysfunction in patients with diffuse glioma: a systematic review of neurocognitive functioning prior to anti-tumor treatment. *J. Neurooncol.* 134, 9–18. doi: 10.1007/s11060-017-2503-z

Varner, V., and Nelson, C. (2014). Cellular and physical mechanisms of branching morphogenesis. *Development* 141, 2750–2759. doi: 10.1242/dev.104794

Wang, Y.-Y., Liu, S.-C., Yang, Z., and Zhang, T. (2010). Impaired hippocampal synaptic plasticity in C6 glioma-bearing rats. *J. Neurooncol.* 103, 469–477. doi: 10.1007/s11060-010-0447-7

Wefel, J., Noll, K., and Scheurer, M. (2016). Neurocognitive functioning and genetic variation in patients with primary brain tumors. *Lancet Oncol.* 17, e97–e108. doi: 10.1016/S1470-2045(15)00380-0

Keywords: dendritic spine, Turing instability, reaction-diffusion model, branching morphogenesis, glioma

Citation: Jia Y, Zhao Q, Yin H, Guo S, Sun M, Yang Z and Zhao X (2021) Reaction–Diffusion Model-Based Research on Formation Mechanism of Neuron Dendritic Spine Patterns. *Front. Neurorobot.* 15:563682. doi: 10.3389/fnbot.2021.563682

Received: 19 May 2020; Accepted: 17 May 2021;

Published: 14 June 2021.

Edited by:

Francesco Rundo, STMicroelectronics, ItalyCopyright © 2021 Jia, Zhao, Yin, Guo, Sun, Yang and Zhao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Xin Zhao, zhaoxin@nankai.edu.cn; Zhuo Yang, zhuoyang@nankai.edu.cn

^{†}These authors have contributed equally to this work