Stochastic simulations of self-organized elastogenesis in the developing lung

In the normal lung, the dominant cable is an elastic “line element” composed of elastin fibers bound to a protein scaffold. The cable line element maintains alveolar geometry by balancing surface forces within the alveolus and changes in lung volume with exercise. Recent work in the postnatal rat lung has suggested that the process of cable development is self-organized in the extracellular matrix. Early in postnatal development, a blanket of tropoelastin (TE) spheres appear in the primitive lung. Within 7 to 10 days, the TE spheres are incorporated into a distributed protein scaffold creating the mature cable line element. To study the process of extracellular assembly, we used cellular automata (CA) simulations. CA simulations demonstrated that the intermediate step of tropoelastin self-aggregation into TE spheres enhanced the efficiency of cable formation more than 5-fold. Similarly, the rate of tropoelastin production had a direct impact on the efficiency of scaffold binding. The binding affinity of the tropoelastin to the protein scaffold, potentially reflecting heritable traits, also had a significant impact on cable development. In contrast, the spatial distribution of TE monomer production, increased Brownian motion and variations in scaffold geometry did not significantly impact simulations of cable development. We conclude that CA simulations are useful in exploring the impact of concentration, geometry, and movement on the fundamental process of elastogenesis.

In the normal lung, the dominant cable is an elastic "line element" composed of elastin fibers bound to a protein scaffold. The cable line element maintains alveolar geometry by balancing surface forces within the alveolus and changes in lung volume with exercise. Recent work in the postnatal rat lung has suggested that the process of cable development is selforganized in the extracellular matrix. Early in postnatal development, a blanket of tropoelastin (TE) spheres appear in the primitive lung. Within 7 to 10 days, the TE spheres are incorporated into a distributed protein scaffold creating the mature cable line element. To study the process of extracellular assembly, we used cellular automata (CA) simulations. CA simulations demonstrated that the intermediate step of tropoelastin self-aggregation into TE spheres enhanced the efficiency of cable formation more than 5-fold. Similarly, the rate of tropoelastin production had a direct impact on the efficiency of scaffold binding. The binding affinity of the tropoelastin to the protein scaffold, potentially reflecting heritable traits, also had a significant impact on cable development. In contrast, the spatial distribution of TE monomer production, increased Brownian motion and variations in scaffold geometry did not significantly impact simulations of cable development. We conclude that CA simulations are useful in exploring the impact of concentration, geometry, and movement on the fundamental process of elastogenesis.

Author summary
The dominant cable in a healthy lung is made up of elastin fibers bound to a protein scaffold and helps maintain alveolar shape and lung volume during exercise. Recent research on postnatal rat lungs indicates that elastogenesis, the development of this cable, is a selforganized process within the extracellular matrix. Elastogenesis is difficult to study by conventional approaches but can be usefully explored by cellular automata simulations. First, elastogenesis involves discrete interactive variables-such as TE monomers and TE spheres-that bind to a protein scaffold. Conceptualizing large mutually interactive populations, particularly difficult in the compressed timeframe of the cable development, is enhanced with simulations. Second, the length scale and geometry of our simulations is a reasonable reflection of the primary septa and the protein scaffold. This scale consistency suggests that the short-range planar interactions in our simulations will faithfully contribute to the large-scale 3-dimensional patterns observed in the developing lung. Finally, elastogenesis appears to be a self-organized process independent of supervised interactions. The simple rules explored in these simulations may not only reflect normal development but also suggest perturbations potentially associated with congenital lung diseases.

Introduction
Elastin fibers are ubiquitous in the extracellular matrix (ECM) of vertebrate tissues. Elastin fibers are an important structural component of the skin, lungs, tendons, cartilage, and cardiovascular system [1]. The elasticity, stability, and durability of elastin fibers are essential for complex functions such as lung ventilation [2]. In the normal lung, the dominant cable is an elastic "line element" composed of elastin fibers bound to a protein scaffold [3]. The cable line element originates in the hilar airways and ends in the subpleural lung [3] (Fig 1). Weibel has described the line element as an "ingenious" fiber continuum that supports the conducting airways as well as the fragile septa of the alveolar walls [4][5][6]. In children, defects in the cable line element can lead to the structural dysfunction associated with bronchopulmonary dysplasia [7]. In adults, disruption of the cable leads to the anatomic and functional limitations associated with emphysema [8].
The cable line element plays a central role in the alveolar phase of lung development. Alveolarization is the process that forms secondary alveolar septa by the lifting new tissue ridges from primitive primary septa. In rats, a species in which alveolar septation occurs postnatally [9], the cable line element forms between postnatal day 4 (P4) and day 14 (P14). Valenzuela and colleagues demonstrated that the process is characterized by a blanket of tropoelastin (TE) spheres detectable in the primary septa on postnatal day 4 (P4) [10] (Fig 2). The ubiquitous spheres had a mean diameter of 2 um and were uniformly distributed in primary alveolar septa [10]. By P14, the tropoelastin appeared to be incorporated into the mature cable line element in a process termed elastogenesis.
An intriguing observation was the spatial relationship of tropoelastin spheres and the protein scaffolding within the primary septa of rat pups [10]. Between P4 and P14, light and electron microscopy demonstrated tropoelastin spheres and elastin monomers linked to the protein scaffold. In contrast, Valenzuela and colleagues found no consistent relationship between the tropoelastin spheres and parenchymal cells. Reminiscent of the observations of elastin self-assembly in vitro [11,12], these findings indicated that the cable was self-organized by physical and chemical processes within the extracellular matrix.
To explore the fundamental process of elastogenesis, we employed cellular automata (CA) simulations. The simulations explored the concentration, geometry and movement of topoelastin spheres in the extracellular assembly of the cable line element.

Elastin aggregation and protein binding
In the early postnatal period (P1-P4), TE monomers in the extracellular space can bind to the protein scaffold either before or after self-aggregation into TE spheres [10]. In our simulation (S1 Movie), a plateau was observed reflecting a dynamic balance between TE monomer generation rate and TE binding to the protein scaffold (both before and after self-aggregation) ( Fig 3A). After the plateau, there was an expected increase in scaffold-associated TE and a commensurate decline in free TE monomer concentration. The average size of the TE spheres decreased as the available TE monomers were bound. The efficiency of TE binding to the protein scaffold suggested an 800-time step process interval; that is, a time interval compatible with the experimental observations of 8-10 days for the development of the mature cable [10]. Importantly, if TE monomers did not aggregate into TE spheres-and the process relied upon TE monomer binding alone-the formation of the cable line element required more than 45 days. This observation supported the role of TE spheres as critical intermediates in the formation of the cable line element.

Spatial distribution of TE monomer production
To assess the impact of the spatial distribution of TE monomer production, our simulation varied the spatial distribution of TE monomer production over a 20-fold range (Fig 4 and S2 Movie). The TE monomer production rate was varied so that the total number of TE monomers was kept constant. The more concentrated production of TE monomers demonstrated a minimal peak preceding a plateau (Fig 4A, black); however, there was little difference between spatial production conditions ( Fig 4B). The spatial distribution of TE monomer production-analogous to the spatial distribution of cells secreting TE monomers-appeared to have only a modest impact on the efficiency of elastin cable formation.

Temporal distribution of TE monomer production
Experimental data indicated that the TE spheres appeared between P1 and P4 [10]. To assess the impact of the rate of TE monomer production, we simulated TE monomer production over a 4-fold range of uniform rates (S3 Movie). Production rates were assessed for their impact on the total area of TE monomer and TE spheres unattached and attached to the protein scaffold ( Fig 5). The highest production rate (shorter production interval) demonstrated an early aggregate plateau and rapid decline ( Fig 5A). In contrast, lower production rates (longer production interval) had a prolonged plateau. Notably, the higher production rate had a supra-additive effect on scaffold binding suggesting the beneficial impact of accelerated production rates (Fig 5C). Qualitatively similar results were obtained with a variable (Gaussian) rate of TE monomer production (S1 and S2 Figs and S4 Movie). The efficiencies gained from higher production rates suggests an adaptive advantage for the burst of elastin transcriptional activity observed in the postnatal period [13].

Brownian motion
The potential impact of self-diffusion, the random motions of TE monomers, was assessed over a 2.5-fold range (S5 Movie). The impact of Brownian motion on the efficiency of TE monomer aggregation ( Fig 6A) and scaffold binding (Fig 6B) was limited. The modest impact was consistent irrespective of the rate of TE monomer production.

TE aggregation affinity
TE monomer aggregation and scaffold binding affinities reflect the physicochemical milieu of elastogenesis. In our simulation, this aggregation and binding affinity was reflected by variation in the critical density of particle aggregation and scaffold binding. A low critical density for aggregation and binding was analogous to high binding affinity. Our simulations were analyzed over a 2-fold range of critical density (S3 and S4 Figs, and S6 Movie). A low critical density (high binding affinity) was associated with a limited number of TE monomersreflecting the rapid increase in TE spheres (Fig 7C). In contrast, high critical density was associated with a high concentration of TE monomers. These results suggest that impaired binding affinity (high critical density) would likely be clinically associated with patchy or incomplete elastogenesis.

Scaffold geometry
The primary septa provide an anatomic boundary for the development of the cable line element. To explore the potential importance of protein scaffold geometry, we tested a variety of scaffold shapes (S5 Fig and S7 Movie). The geometric shapes had no effect on the efficiency nor the appearance of the final scaffold. These results suggest that our findings are not limited by our idealized scaffold geometry.

Discussion and conclusion
Previous work has shown that cable line element is assembled in the extracellular spaceapparently independent of cellular instruction [10]. Here, we use cellular automata simulations to show that the extracellular assembly of the cable line element is consistent with a self-organized process dependent upon the spatial and temporal distribution of extracellular TE spheres. Elastogenesis is a fundamental developmental process that is difficult to study by conventional approaches but can be usefully explored by cellular automata simulations. First, elastogenesis involves discrete interactive variables-such as TE monomers and TE spheres-that bind to a protein scaffold. Conceptualizing large mutually interactive populations, particularly difficult in the compressed timeframe of the cable development, is enhanced with simulations. Second, the length scale and geometry of our simulations is a reasonable reflection of the primary septa and the protein scaffold. This scale consistency suggests that the short-range planar interactions in our simulations will faithfully contribute to the large-scale 3-dimensional patterns observed in the developing lung. Finally, elastogenesis appears to be a self-organized process independent of supervised interactions [10]. The simple rules explored in these simulations may not only reflect normal development [14,15] but also suggest perturbations potentially associated with congenital lung diseases [16]. An intriguing observation from these simulations was the impact of the TE spheres on the efficiency of cable development. Secreted from elastogenic cells as a 60-70 kDa protein, TE is the soluble precursor of elastin. The TE monomer has alternating hydrophobic and hydrophilic domains [17]. The interactions between hydrophilic and hydrophobic domains-a process called coacervation-facilitates the self-assembly of TE into quantized spheres [18]. The TE spheres have been postulated to function as intermediates in elastin macroassembly [19]. Lysine residues on the TE sphere surface facilitate the lysyl oxidase-dependent cross-linking of TE spheres to not only to other spheres, but also extracellular scaffold proteins [20]. Our simulations suggested an adaptive advantage of the process of coacervation and the development TE sphere intermediates; namely, the formation of TE spheres enhanced the efficiency of elastin binding to the protein scaffold by an estimated 5-fold.
A limitation of this work is our simplification of the experimental process of structural selforganization. The initial steps (1 and 2) in the process-the generation of TE monomers and aggregation of TE monomers-are below the limits of our microscopic resolution. Nonetheless, we know from molecular data that TE monomer transcription is active during this phase [21]. A second limitation is our simplified assembly of the cable line element. The protein scaffold involves numerous extracellular components. A growing list of more than 30 different molecules, identified by microscopy and immunohistochemical approaches, have been found to interact with elastic fibers and microfibrils [22]. The interactions of these extracellular components is an important frontier in understanding the process of self-assembly.
The random and uniform distribution of the TE spheres at the onset of secondary septation suggests that it is the scaffold proteins, not the tropoelastin spheres, that are responsive to mechanical force distribution [10]. Tomoda et al. have shown that the collagen fiber orientation closely correlates with respiratory movement [23]. Pins and colleagues as well as Chow and coworkers have shown that the axial alignment of collagen fibers can be achieved by mechanical loading of the uncrosslinked fibers [24][25][26]. Fibrillin microfibrils also undergo gross changes in molecular configuration when subjected to stretch [27]. Although less wellstudied than collagen, fibrillin microfibrils subjected to stretch undergo significant changes in molecular conformation [27]. The conformational effects of mechanical loading provides a useful explanation for the organ-sized length scales of the cable line element; that is, extracellular loads contribute to the efficiency of a structure that extends from the central airways to the subpleural matrix [3].
The results of our simulations-based on in vivo observations [10]-suggests that the cable line element self-organizes in the extracellular matrix. Although elastogenic cells such as fibroblasts, muscle cells and chondrocytes, produce tropoelastin, the physicochemical properties of the extracellular matrix contribute to tropoelastin formation and scaffolding interactions. The observations are reminiscent of in vitro elastin self-assembly [11,28]. We speculate that lung ventilation provides a motive force for this process. Finally, self-organized elastogenesis has practical implications for normal lung development and the pathologic processes of bronchopulmonary dysplasia. We anticipate this will be a promising area of future research.

Cellular automata simulation procedures
According to the experimental observations (Figs 1 and 2), the process of the structural selfassembly of collagen and elastin in the developing mammalian lung can be divided into the following steps: generation of the TE monomer; TE monomer aggregating into TE spheres; aggregation of TE spheres; Attachment of the TE spheres to the protein scaffold. We use cellular automata to simulate these steps. The code is implemented using MATLAB (2021a). The detailed simulation procedures are explained below.

Generation of tropoelastin monomer
TE spheres are formed by the self-assembly of TE monomers. The TE monomers are generated in a limited time (denoted by generate_Time) in the simulation. As the size of TE monomers is around 15 nm, we used a value of each square unit in the simulation domain to represent the contents of the tropoelastin monomer. We use red color to represent the TE monomers and use transparency to represent the value magnitude. A lower transparency represents a higher value. In each step within gener-ate_Time, the value of each unit increases by generate_Tol with the probability of generate_rate. Two types of generation values were used: a constant generation value or a Gaussian distributed generation value. The constant or Gaussian "generation values" distributions are over time. The total generated tropoelastin monomers are the same in the simulations with Gaussian distributions and with constant generate value.
As the TE monomers and TE undergo random Brownian motions, Brownian motions are included in the simulation. The Margolus neighborhood was used to mimic the Brownian motion of the TE monomer. The particle movement strategy is shown in S6 Fig

TE monomoer self-assembly into TE spheres
The tropoelastin monomer will undergo self-assembly to form tropoelastin. To model this process, the tropoelastin monomers are transited to tropoelastin spheres under two situations: (1) the value representing the monomer's content reaches a threshold σ 1 , or (2) the values of the Moore neighborhood (composed of a central cell and the eight cells that surround it) added together to reach a threshold σ 2. Both situations are based on the reasonable hypothesis that tropoelastin spheres are formed when the content of the monomers exceeds a nucleation threshold [29]. In the simulation, each light blue sphere represents a TE. As the minimum diameter of the tropoelastin sphere is around 1 μm, the length of one unit is set as 1 μm. Brownian motion is also applied to the tropoelastin spheres.

Aggregation of tropoelastin
From the experimental observation that the diameters of tropoelastin spheres span from 1-6 μm (Fig 2), we can conclude that the tropoelastin spheres will aggregate to form larger coalesced tropoelastin spheres. In the simulation, if two tropoelastin spheres contact, they will adhere and form larger tropoelastin spheres. A division principle is used to prevent the overgrowth of the tropoelastin. If the size of a tropoelastin sphere is larger than δ 1 , it will split into two spheres, similar to the division of large water droplets. In the simulation, we choose δ 1 = 6.5 2 = 42, corresponding to a full-sized coacervated TE with the diameter of around 6-7 μm, which agree well with the experimental observations that the elastin of the postnatal rats was detected as 1-6 μm tropoelastin spheres. δ 2 represents the minimum size of the TE that can crosslink with the collagen scaffolds. We set δ 2 =~0.9δ 1 = 38, as most crosslinked TE are large from experimental observation.

Attachment of TE to the protein scaffold
The in vivo experiments demonstrate apparent crosslinking of the TE spheres to the protein scaffold (Fig 2). To simulate the adhering process, the following rule is used. If a coalesced tropoelastin is adjacent to the fibrillar structures and its size is larger than δ 2 , it will crosslink to the collagen scaffolds. The pseudocode for the cellular automata simulation are shown in Table 1. The standard parameters used in the simulation are shown in Table 2. Three different cable line shapes are used in the simulation. The cable lines are assumed to be straight lines which forms square, hexagon and triangular shapes, respectively. The linewidth is 4 pixels. The cable lines form square, hexagon or triangular periodic units with the same areas. Periodic boundary conditions are used.

Time and length scale
The time and length scale of the simulation is discussed. As mentioned before, the size of a TE sphere is around 1 μm and is represented by 1 pixel (a unit cell). Therefore, the length scale is λ l = 1 μm/pixel. The length scale in experiments and simulations is shown in Fig 2. the TEM images of the TE of rats on postnatal days 4 are shown in Fig 2. Azure blue staining provides sufficient contrast to demonstrate the broad distribution of the sphere. Fig 2F shows the simulation image of TE. The computation domain is chosen as 200×200 μm 2 . It can be seen that the length scale of the simulation agrees with the experimental observations. Next, the time scale is studied. For a random walk, the mean squared displacement of a particle is proportional to the time interval [30]: where D is the diffusion coefficient, λ t is the time scale. The diffusion coefficient of a particle can be calculated using the Stokes-Einstein equation in Brownian motion where k B is the Boltzmann's constant. T is the absolute temperature, η = 1.9 Pa�s is the viscosity coefficient of the surrounding liquid [31]. d = λ 1 is the particle diameter. According to the above Stokes-Einstein equation, small particle size, low viscosity of the surrounding fluid and high temperature result in faster motion. D is calculated as 2.412×10 −16 m 2 /s. The time scale can then be calculated using Eq (1) as λ 1 �2072 seconds = 0.024 days. Experiments show that a blanket of tropoelastin (TE) spheres is detectable in the primary septa on postnatal day 4 (P4) (Fig 2), and the tropoelastin appeared to be incorporated into the mature cable line element by P14. Therefore, the generate_Time of the tropoelastin monomers is between days 4 and 14.
We choose the generate_Time as 500 in the simulation, which corresponds to 12 days. Brownian motion simulations of multiple particles are conducted, as shown in S7 Fig