Quantitative modeling of regular retinal microglia distribution

Microglia are resident immune cells in the central nervous system, showing a regular distribution. Advancing microscopy and image processing techniques have contributed to elucidating microglia’s morphology, dynamics, and distribution. However, the mechanism underlying the regular distribution of microglia remains to be elucidated. First, we quantitatively confirmed the regularity of the distribution pattern of microglial soma in the retina. Second, we formulated a mathematical model that includes factors that may influence regular distribution. Next, we experimentally quantified the model parameters (cell movement, process formation, and ATP dynamics). The resulting model simulation from the measured parameters showed that direct cell–cell contact is most important in generating regular cell spacing. Finally, we tried to specify the molecular pathway responsible for the repulsion between neighboring microglia.


Model of microglia process length distribution
Based on the experimentally observed motion of the process, we modeled the dynamics of process length distribution u(x,t) using advection and degradation terms as follows: g is the velocity of process extension and p is the collapse probability of process. By the model (5), the steady-state of length distribution u(x, •) can be described as follows: U 0 represents a certain constant. u(0, •) is exponential distribution, indicating there is a linear relationship between the log 435 of process number and process length.

436
The sum of the process length is described by the definite integral as follows: U 0 are obtained by linear fitting equation (U 0 = 7.6 in avascular area and U 0 = 7.7 in vascular areas). U 0 g 2 p 2 are 1.0 ⇥ 10 5 µm 437 vs. 1.1 ⇥ 10 5 µm, and the sum of the process length from direct measurements are 1.1 ⇥ 10 5 µm vs. 1.3 ⇥ 10 5 µm, respectively. 438 This similarity shows that the model captures the characteristics of the actual process distribution.

440
• c: We calculated the ATP chemotaxis coefficient using the published data of ATP chemotaxis of microglia in the Dunn 441 chamber 8 .

442
• h: We determined the ATP decay coefficient by fitting the data that ATP decay of homogenized retina 32 . The fitting equation was as follows: • p: The basal extracellular ATP production rate Sup1 is calculated by steady state concentration of ATP (a 0 ) and ATP decay 443 rate h. We used extracellular ATP of the bovine retina (6.2 ± 0.7 nM Sup2 ) as a 0 . We obtain p by steady state equilibrium 444 0 = p ha 0 . Uptake by microglia is small and neglected in this calculation. We measured the average NND (31 ± 0.12 µm, n = 7053, Fig. 3) in P5 mouse retina; thus, we assumed microglia 447 repulsion radius is 45 µm, considering the ratio of their report 17 .

448
• a w : Intracellular ATP is millimolar order 33 . We set production of ATP at the wound site so that a w becomes 1/300 of 449 intracellular ATP.

450
Screening repulsion-related genes from single-cell RNASeq data. 451 It is known that there is a regional difference in microglia cell density. Microglial cell density is lowest in the cerebellum and 452 highest in the hippocampus 57 . This may be correlated to the mutual repulsion of microglial cells.

453
At first, we obtained differentially expressed genes of microglia from different regions of CNS (GSE123025) 60 . DEGs 454 between the cerebellum and other regions of CNS were provided in supplemental files 60 (FDR < 0.05).

455
• Genes upregulated in cerebellum: Again, gene enrichment analysis using these DEGs did not show cell movement-related pathways. However, among these 463 DEGs, Egr1 and Klf2 were detected in both cases. They are relatively well studied and are not general housekeeping genes. 464 This information may facilitate further candidate-approach studies.