Remodeling mechanisms determine size distributions in developing retinal vasculature

The development of retinal blood vessels has extensively been used as a model to study vascular pattern formation. To date, various quantitative measurements, such as size distribution have been performed, but the relationship between pattern formation mechanisms and these measurements remains unclear. In the present study, we first focus on the islands (small regions subdivided by the capillary network). We quantitatively measured the island size distribution in the retinal vascular network and found that it tended to exhibit an exponential distribution. We were able to recapitulate this distribution pattern in a theoretical model by implementing the stochastic disappearance of vessel segments around arteries could reproduce the observed exponential distribution of islands. Second, we observed that the diameter distribution of the retinal artery segment obeyed a power law. We theoretically showed that an equal bifurcation branch pattern and Murray’s law could reproduce this pattern. This study demonstrates the utility of examining size distribution for understanding the mechanisms of vascular pattern formation.


Introduction
The development of the retinal vasculature has extensively been investigated as a model system for vascular pattern formation [1]. In mice, the retinal vasculature is established during the perinatal period. Initially, astroglia migrates from the optic disc region, developing an astroglial meshwork. Next, the endothelial cells migrate from the optic disc region on this preexisting astroglial meshwork, forming the first capillary network beginning at postnatal day 1 (P1) (Fig 1a). Subsequently, the hyaloid arteries connect to multiple sites within the capillary network near the optic disc, thereby dramatically increasing the blood flow to the network and driving vascular remodeling. Initially, blood with high oxygen concentration runs through the arteries, leading to capillary regression around the arteries. This region of capillary retardation is called the avascular zone [2]. Simultaneously, the blood flow causes changes in the arterial segment diameter, resulting in the formation of an arterial vascular tree [3]. Development At first, the segment radii of the capillary network was random at P1 (Fig 1b). After the vascular flow was established, the arterial tree structure radiating from optic disc became evident at P4 (Fig 1c). Finally, the arterial trees became distinguishable at P7 (Fig 1d, red arrows). Various theoretical models have been proposed to understand the vascular pattern formation. Classically, the general relationship between the vascular radius before and after branching has been explained in terms of an energy minimization problem known as Murray's law [4]. Similarly, the general argument of scaling in various networks has been proposed, which has become an influential concept in physics [5]. More recent studies have reproduced the retinal vasculature remodeling process observed in the avian yolk sac [6] and retinal vasculature [7] using a simple vascular diameter growth rule. Additionally, a large scale computational model that combined the molecular regulation of cell migration, vascular flow, and gas exchange succeeded in reproducing the pattern formation of retinal vasculature [8].
Quantification of the retinal vascular structure has been done extensively for understanding mechanisms of pattern formation. Previous efforts concentrated on quantifying various aspects of the vascular pattern [9,10]. Using image processing techniques, various quantities such as vascular outgrowth, vascular area fraction, vessel segment length, number of branches, etc. were measured at a sequential stage of development to obtain characteristics of the vascular network shape [10].
However, the relationship between theoreitical models and measured quantities has not yet been thoroughly addressed. For example, the fractal dimension was frequently used as a measure of the vascular pattern [11,12]. In these studies, the fractal dimension of retinal vasculature was quantitatively measured using the mass-radius or box-count method. However, it is unclear how these quantities are correlated to the pattern formation mechanism. A previous study by [11] asserted that the pattern formation mechanism is based on either diffusion-limited aggregation or Laplacian growth owing to the similarity of the fractal dimension D = 1.7. However, based on our current knowledge of the above described retinal vasculature development, this argument is too naive.
In the present study, we combined the above mentioned recent experimental findings, mathematical modeling, and image processing techniques to elucidate the relationship between vasculature size distribution laws and pattern formation mechanisms. We focused on two pattern formation processes-(1) formation of the avascular zone around arteries, and (2) vascular remodeling of arterial trees. First, we focused on the formation of the avascular zone, and experimentally observed that the island size distribution obeyed an exponential law. We formulated a minimal model for the stochastic disappearance of capillaries surrounding the arteries, and theoretically demonstrated that the island size distribution generated by the model exibited the same exponential distribution as was observed experimentally. Second, we focused on arterial remodeling, and experimentally found a power distribution of the radii of the arterial vascular tree. We also confirmed that Murray's law was established during the development of retinal vasculature. Thereafter we theoretically showed that the combination of bifurcated geometry and Murray's law could result in power distribution between the segment radius and the number of arterial vascular segments. These results show the utility of examining size distribution for understanding the mechanisms of pattern formation in retinal vasculature.

Measurement of the island size distribution
We captured the images of P1 to P8 mouse retinas. For image processing, we used the Fiji software [13]. First, we sharpened the image, transposed it to an 8-bit grayscale image, and manually determined the threshold. Subsequently, we measured the island size using the "particle analysis" command. The distribution was visualized and analyzed using Mathematica (Wolfram Research Inc.).

Measurement of vessel diameter using distance map
We measured vessel diameters by combining the distance map and the skeletonized image (Fig 4). First, we made a binarized image of the arterial region using aSMA and PECAM images by "AND" operation, background subtraction, and thresholding (a). After that, we skeletonized the binarized image (b). Subsequently, we detected the three-way junction points of this skeletonized image using 3 × 3 averaging kernel (c). Removal of these junction points (c) from the skeletonized images (b) allowed us to generate a skeletonized image of vessel segments (e). Following this, we created a distance map using the same binarized image (d). Finally, we multiplied the skeletonized image of vessel segments and the distance map and measured vessel segment diameters using the "Particle analysis" command. Short segments (<3 pixels) and extremely thick segments (>7 pixels) were ignored. Image processing was performed using Fiji software [13], whereas quantification was performed using Mathematica. Additionally, we assessed the effect of lattice on segment rotation, and the error was <5% (S3 Fig in S1 File).

Island size distribution in the vascular network
Quantification of the island size distribution in the retinal vasculature. First, we measured the island size distribution in the retinal vasculature (Fig 2). Vascular structures were detected using the IB4 lectin stain, and the positively stained region was binarized by thresholding. Subsequently, the arterial regions were skeletonized, and the island size distribution was automatically detected using the Fiji software (Fig 2b and 2c). The log-linear plot of the histogram showed linearity in size distribution (Fig 2d. A distribution of one representative sample is shown for each time point.). All data points fit well to the linear model, with a exceptions for the smallest island size group in each plot (Fig 2d, red circles). Moreover, we measured the time course of island size distribution and found that observed linearity was established as the development of the avascular zone (Fig 2d).
A model to generate large islands in the avascular zone surrounding the arteries. In this section, we propose a simple model that produces islands that exibit exponential size distribution. It has been known that with the establishment of arterial flow, the capillaries near the arteries disappear, owing to the increased oxygen concentration (Fig 3a, [2]). As a result, llarge islands are distributed near arteries close to the optic disc. The region region called the avascular zone. To model this phenomenon, we assumed that the initial vascular meshwork was a regular rectangular lattice, and implemented the well-known fact that the capillaries directly connecting to the arteries disappear in a stochastic manner. In the model, the island distribution is established by the following two-step procedure: 1. A regular mesh of a capillary network is generated.
2. The inter-island capillaries adjacent to the arteries disappear by the probability p due to high oxygen concentration.
For simplicity, we assume the initial capillary network is a rectangular lattice. Using this model, we analytically obtained the frequency of islands of size n after the stochastic disappearance of capillaries (Table 1). The lattices which are not adjacent to the arteries (N far ) do not change, resulting in the large number of small islands as observed in Fig 2d (red circles).
In this model, the frequency of islands of size n was roughly proportional to p n−1 , resulting in a linear pattern in the log-linear plot. We also confirmed this tendency numerically (Fig 3c). Therefore, this model explains the observed exponential distribution of the island size in the retinal vasculature.

Power distribution of vascular diameters in the retinal artery network
Diameter distribution measurement in the retinal vasculature. Next, we sought to experimentally assess whether the vasculature exhibited Murray's law and a power distribution between thickness and frequency of vessel segments in the developing mouse retinal arteries. Segment radii were automatically measured using Mathematica (Fig 4). Stage First, we chose appropriate stages for retinal artery differentiation (P5-8) and obtained arterial diameter using aSMA and PECAM stain (Fig 4). To confirm whether Murray's law holds in the retinal arterial trees, we plotted r 3 1 þ r 3 2 and r 3 0 , with r 0 indicating the thickest segment. At developmental stages we observed, we obtain a strong correlation of the two variables in the plot (Fig 5b, R 2 > 0.9 in all stages), indicating Murray's law was established by the observed time points. In addition, we plotted the frequency and the thickness of vessel segments in a log-log plot. The power distribution between thickness and frequency was established by stages P7 and P8, as indicated by the linearity of the graph (Fig 5c).

PLOS ONE
We additionally examined whether the branch angles exhibited Murray's law, as has previously been predicted [14]. However, the experimentally measured angles did not fit the law (S2 Fig in S1 File).
Relationship between Murray's law and vascular branch power distribution. In this section, we show that the combination of Murray's law and stochastic bifurcation of the retinal artery can result in a power distribution between a vessel segment number and diameter [5,14,15]. In Murray's law ( [14], Fig 6), the following relationship holds at the three-way  Table 1. Expected values of number of clusters. N near , N far and p represents the initial number of lattices around arteries, the initial number of lattices that are not adjacent to the arteries, and extinction probability of vasculature between islands near arteries respectively. junctions of blood vessel:

Cluster size Number
Assuming that the retinal artery tree consists of repeated equal bifurcations, we obtain two proportional relations-one between the mother vessel radius r n and the sister vessel radius r n+1 , and the other between the number of mother vessels N n and that of sister vessels N n+1 . Therefore, the number and radius should obey a power law. With the progression to the next generation, the number of n generations doubles from that of n − 1th generation, and the radius becomes 2 −1/3 times that of the former generation (Fig 6b). Therefore, when drawing a log-log plot, the gradient becomes −3 (Fig 6c). The establishment of Murray's law can be explained by setting an appropriate vessel diameter growth function (S1 Text, S1 Fig in S1  File).
With a probabilistic equal bifurcation, the log-log plot exhibits a gentle gradient (Fig 6d). When q is the probability of generating the bifurcation, we can estimate the gradient of the plot accompanying a change of q as follows:  : [14]. (b) When the tree structure is generated by equal bifurcation under Murray's law, it is possible to explicitly obtain the number and radius of the n + 1th bifurcation. (c) Each blood vessel segment radius and number should show linearity on a log-log plot, and the gradient should be −3 [16]. (d) Incomplete bifurcation. In the real retina, the arterial tree did not form a complete bifurcation tree. If we assume that the bifurcation occurs at probability q, the number of thin branches decreases, resulting in a smaller gradient. (e) Relationship between bifurcation probability q and gradient. When q is less than 1, the gradient is larger than −3.

Discussion
In the present study, we evaluated the relationship between size distribution laws and pattern formation mechanisms. Among the size distribution laws, power distribution could be correlated to the fractal dimension. Following the original study by [11], several similar investigations were performed focusing on a fractal dimension of retinal vasculature as a diagnostic tool [12,[17][18][19][20][21][22]. To date, only a small number of studies have correlated size distribution with fractal pattern formation mechanisms, which were not based on the current understanding of retinal vasculature development [23].
Since the two models we presented reflect the specific biological processes ((1) disappearance of vessel segments by high oxygen concentration and (2) arterial remodeling by flow), experimentally obtained p and q can be used as a measure for defective pattern formation in the retinal vasculature. For example, p can be correlated to the degree of oxygen supply during retina development. q can be interpreted as a measure for vascular remodeling by flow and arterial differentiation. It would be intriguing to use these parameters to investigate the shape change induced by various factors such as oxygen-induced retinopathy or absence of pericytes [24]. Strain difference of retina vasculature has been reported [10], and it is intriguing to detect the difference using these parameters since they are directly correlated to pattern formation mechanisms.
A previous report has confirmed the prediction of Murray's law by manually measuring pig coronary arteries [9,25]. A modification of Murray's law has been proposed by [26], which considers the vessel wall metabolic cost as follows: In this equation, K w represents the metabolic constant, and h indicates the vessel wall thickness. In this case, when dE dr ¼ 0 we obtain f / ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi r 6 þ cr 5 p ð6Þ where c / h. Thickness Since we observed the wall thickness existed at this stage (S4b Fig in S1 File), this effect might influence the result. The biological mechanisms of endothelial regression by high oxygen concentration and of vessel diameter change by flow have been well studied. High oxygen concentration results in the inhibition of vascular endothelial growth factor (VEGF) via hypoxia-inducible factor (HIF) [27]. As a result, capillaries are not maintained due to the lack of VEGF. In addition, it has been shown that vessel diameters are influenced by blood flow [28] and pericyte activity [29,30]. Recent studies have suggested that shear stress regulates endothelial and smooth muscle cell signaling. This signaling is mediated by several factors, including nitric oxide (NO), prostaglandin I-2 (PGI-2), platelet-derived growth factor (PDGF-BB), transforming growth factor -β1 (TGF-β1), and microRNAs (miRs) [31]. There are both fast (NO and PGI-2) and slow (PDGF-BB, TGF-β1, and miR126) factors, which are used for communication between the endothelial and smooth muscle cells [32,33]. Owing to our interest in this structural change, we focused on the latter factors. PDGF-BB and TGF-β1 are produced by endothelial cells under low shear stress. Functionally, PDGF-BB activates smooth muscle proliferation, migration, and contraction, whereas TGF-β1 activates smooth muscle cell differentiation. Additionally, miR126, which is produced by endothelial cells under laminar shear stress, activates smooth muscle cell proliferation. In future studies, these slow factors can be used to assess the relationship between flow and diameters experimentally. Pericytes are also known to be involved in the retinal vasculature remodeling process. It has been shown that pericytes play a role in modulating the extracellular matrix degradation, production, and assembly via an interaction with endothelial cells [29,30]. Goal Manipulation of pattern formation and size distribution using these molecules, or application of a similar method to other vascular systems like that in the central nervous system would be an interesting future goal.