On the preservation of vessel bifurcations during flow-mediated angiogenic remodelling

During developmental angiogenesis, endothelial cells respond to shear stress by migrating and remodelling the initially hyperbranched plexus, removing certain vessels whilst maintaining others. The key regulator of vessel preservation is cell decision behaviour at bifurcations. At flow-convergent bifurcations where migration paths diverge, cells must finely tune migration along both possible paths if the bifurcation is to persist. Experiments have demonstrated that disrupting the cells’ ability to sense shear or junction forces transmitted between cells impacts the preservation of bifurcations during the remodelling process. However, how these migratory cues integrate during cell decision making remains poorly understood. Therefore, we present the first agent-based model of endothelial cell flow-mediated migration suitable for interrogating the mechanisms behind bifurcation stability. The model simulates flow in a bifurcated vessel network composed of agents representing endothelial cells arranged into a lumen which migrate against flow. Upon approaching a bifurcation where more than one migration path exists, agents refer to a stochastic bifurcation rule which models the decision cells make as a combination of flow-based and collective-based migratory cues. With this rule, cells favour branches with relatively larger shear stress or cell number. We found that cells must integrate both cues nearly equally to maximise bifurcation stability. In simulations with stable bifurcations, we found competitive oscillations between flow and collective cues, and simulations that lost the bifurcation were unable to maintain these oscillations. The competition between these two cues is haemodynamic in origin, and demonstrates that a natural defence against bifurcation loss during remodelling exists: as vessel lumens narrow due to cell efflux, resistance to flow and shear stress increases, attracting new cells to enter and rescue the vessel from regression. Our work provides theoretical insight into the role of junction force transmission has in stabilising vasculature during remodelling and as an emergent mechanism to avoid functional shunting. Author Summary When new blood vessels are created, the endothelial cells that make up these vessels migrate and rearrange in response to blood flow to remodel and optimise the vessel network. An essential part of this process is maintaining the branched structure of the network; however, it is unclear what cues cells consider at regions were vessels branch (i.e., bifurcations). In this research, we present a computer model of cell migration to interrogate the process of preserving bifurcations during remodelling. In this model, cells at bifurcations are influenced by both flow and force transmitted from neighbouring cells. We found that both cues (flow-based and collective-based) must be considered equally in order to preserve branching in the vessel network. In simulations with stable bifurcations, we demonstrated that these cues oscillate: a strong signal in one was accompanied by a weak signal in the other. Furthermore, we found that these cues naturally compete with each other due to the coupling between blood flow and the size of the blood vessels, i.e. larger vessels with more cells produce less flow signals and vice versa. Our research provides insight into how forces transmitted between neighbouring cells stabilises and preserves branching during remodelling, as well as implicates the disruption of this force transmission as a potential mechanism when remodelling goes wrong as in the case of vascular malformation.


Introduction
Angiogenesis occurs as two distinct phases: an early phase in which sprouting neovessels assemble to form the initial immature vascular plexus, and a late phase which remodels the plexus into its final functional form [1,2]. Blood flow is typically shunted away from sprouting endothelial cells (ECs), which exhibit a dynamic exploratory phenotype as they invade the avascular space. However, ECs exhibit a fundamentally different phenotype upon receiving blood flow: in response to experiencing the physical force of flow (i.e., wall shear stress along the luminal surface, WSS), ECs will re-align their polarity and migrate against the direction of flow [3][4][5][6][7][8]. This phenomenon, referred to as EC flow-migration coupling [8], has been shown to be the primary driver of vascular remodelling during developmental angiogenesis.
There are two principal outcomes of flow-mediated remodelling: the pruning of superfluous connections that arise during the sprouting phase, and the establishment of vascular hierarchy via diameter control. Under the flow-migration hypothesis, ECs move from lowflow vessels as they are attracted to high-flow vessels along a path opposite to flow. This attraction results in pruning of inefficient vessels while reinforcing the established preferred flow paths, including precursors to arteries and veins. This hypothesis presents a paradigm shift through which we can now view quiescent vascular structures as nonlinear dynamic systems of migrating cells which have reached stability at a critical point and shifts between healthy vascular tissue and diseased now as transitions between stable fixed points which may be reversible. However, the fundamental question as to how vessel networks determine the appropriate structure and hierarchy needed to establish optimal tissue perfusion during flowregulated migration and remodelling remains largely unanswered.
Branching within the vasculature is a key component to effective transport, and it appears that an important result of flow-mediated remodelling is to remove some complexity from the initial plexus in order to improve efficiency while maintaining enough complexity within the network to ensure large surface area and short diffusion distance to the surrounding tissue. However, under the basic flow-migration hypothesis, in which ECs simply move against flow from areas of low-flow to high, a discrepancy arises at bifurcations: what do ECs do at flow-convergent bifurcations in which two paths to migrate against flow exist? Do ECs simply choose the high-flow branch, thereby reinforcing these vessels at the expense of their low-flow counterparts? Or does some mechanism exist that allows bifurcations with a flow difference between the two branches to persevere during remodelling? The flow-migration hypothesis is relatively new in its development and lacks a description of EC behaviour at vessel bifurcations as well as consideration as to how certain bifurcations are preserved and not others.
Very little is known about the mechanisms regulating EC flow-migration coupling; however, some signalling pathways have been revealed to be involved. In particular, the noncanonical Wnt signalling pathway plays a profound role. Franco et al. demonstrated that knocking down Wnt5a in mutant mice (hence referred to as Wnt KD) rendered ECs in the developing retina more sensitive to WSS, increasing their axial polarisation against flow and migration levels [6]. Additionally, they found reduced amounts of bifurcations and increased regression events within these mice, indicating that increased levels of flow-mediated remodelling reduced the branching complexity of the emerging vasculature. Noncanonical Wnt signalling has also been implicated in sprouting ECs as well in a different role. Carvalho et al. demonstrated the role of Wnt5a in collective cell behaviour at the sprouting front, where Wnt5a works to reinforce leader-follower collective polarity by strengthening and stabilising adherens junctions, facilitating force transmission between ECs [9].
Putting these findings together demonstrates that interfering with noncanonical Wnt5a signalling reduces junction force transmission between ECs, which facilitates better flow-based polarisation and migration, increases vascular remodelling, and results in a network that is less branched. The fact that reducing junction force transmission facilitates flow-based remodelling rather than inhibits may seem counterintuitive, as many collective migration processes depend on junction force transmission [10,11]. Rather, it would seem that flow-migration coupling during remodelling is a different type of process all together, and the "collectiveness" of the ECs interferes with the process in order to keep remodelling in check. This implies that flowmigration coupling is more of an "individualistic" process in which isolated ECs responding to a global signal (i.e., WSS due to blood flow), rather than a process regulated via cell-cell communication.
We hypothesise that EC migration during angiogenic remodelling results from integrating both flow-directed migration (guided by differences in luminal shear stress) and collective migration (guided by cell-to-cell junctional force transmission). Furthermore, we propose that these two migratory cues interact (either cooperatively or competitively) to determine bifurcation stability and promote or avoid vessel regression during remodelling. In this work we present an agent-based model (ABM) of EC flow-migration coupling to demonstrate that competition between flow-based and collective-based cues can stabilise vessel bifurcations during remodelling. ABMs have proven to be useful tools in the study of sprouting angiogenesis due to their naturally discrete representation of cells and their ability to capture emergent behaviour [12][13][14][15][16][17][18][19][20]; however, there are currently no models of the coupling between flow and migration during remodelling to date. In ABMs, we can prescribe various "rules" of ECs behaviour at bifurcations and observe the outcome as an emergent property, allowing us to characterise the impact of different mechanisms when it comes to achieving bifurcation stability. We found that bifurcation stability cannot be achieved when ECs follow flow or collective migration cues alone. The best outcome for bifurcation stability was when we included near equal contributions of shear stress and junction force transmission cues when determining EC migration at bifurcations. Additionally, we found that this stability arises due to the competitive interplay of the two mechanisms, which oscillate back in forth to keep the bifurcation stable. The competitive nature of these two cues is of hydrodynamic origin, arising from the inverse relationship between lumen diameter and vessel resistance to flow and can be characterised by a single parameter. We postulate that published observations on the molecular regulation of EC flow-migration coupling can be re-envisioned under this new lens, offering a powerful rationale for developing a mechanistic understanding of vascular remodelling dynamics at the level of the whole vascular plexus.

Agent-based model of EC migration coupled to flow
In this study we present an agent-based model of EC migration coupled to blood flow within an idealised bifurcated vessel network (hence referred to as the A branch, see the Methods Section for a complete description). This network consists of a feeding vessel and a draining vessel, connected by a proximal branch and a distal branch ( Fig. 1 A). Each vessel was discretised into segments and each segment was seeded with an initial number of "agents" representing ECs (Fig. 1 B). Flow was driven by the difference in pressure prescribed at the inlet and outlet in order to recapitulate previously reported WSS values [21], and calculated via the Hagen-Poiseuille equation while treating blood as a Newtonian Fluid with constant viscosity. The flow conductance (i.e., the inverse of the flow resistance) was calculated using a three-dimensional (3D) approximation of the lumen diameter by wrapping the number of cells in each vessel segment into the circumference of a circle (Fig 1 C). Each simulation was run for a prescribed number of steps in time, and during each step ECs moved against the direction of flow to the neighbouring upstream segment. Periodic boundary conditions were prescribed to handle the case of ECs at the inlet, in which case they migrated to the outlet (thus preserving the total number of cells during each simulation). After each time step, vessel diameter was updated depending on the new number of ECs in each segment and flow recalculated, thus directly linking flow with migration.

Wall shear stress alone cannot stabilise bifurcations
Our initial investigation involved implementing a series of rules for behaviour at the bifurcation (or bifurcations rules, BRs) and observed the emergent outcome in the resulting simulations. In the first bifurcation rule (BR 1), we programmed that upon reaching the flowconvergent bifurcation, the EC would choose the branch with larger shear stress (Fig 2 A, S1 File). This simulation always resulted in the loss of the low-flow distal branch, as cells would always turn into the high-flow branch, meaning cells migrating out of the distal branch were not replenished with new cells. This resulted in a loss of the bifurcation within the network as the cells formed into a single path from inlet to outlet along the proximal branch. Turning into the high-flow branch as in BR 1 means that cells would have to make a sharp change in their migration direction. Therefore, in BR 2 we interrogated if cells would prefer the path that required the smallest change in direction (Fig 2 B, S2 File). This simulation always results in the loss of the proximal branch and reinforcement of the distal branch, as this path requires no direction change from ECs at the bifurcation.
We then asked the question of how the flow-convergent bifurcation may remain in such a network configuration. It would seem that for the bifurcation to remain stable, some of the cells would have to choose the high-flow option, while others choose the path requiring the least change in direction/polarisation. Therefore, we implemented a new mechanism in which each branch was assigned equal probability of drawing in new cells ( , BR 3). In 1 2 0.5 P P   these simulations, both the proximal and distal branch remained stable, preserving the bifurcation, and stabilised at a similar diameter despite the difference in flow between them  the distal path is twice the length of the proximal path, and it takes several trips around before the smoothing algorithm settles down the diameter fluctuations. (C) Mean diameter and standard deviation for the 10 runs using BR3.. Using a random number generator to determine which branch each cell entered with equal probability resulted in stabilisation of both branches and no discernible difference in diameter between the two branches. (D) Mean diameter and results of each run using BR 4. Using unequal probability between the two branches, while favouring the high-flow branch, resulted in a form of diameter control, as the high-flow proximal branch stabilised at a larger diameter than the low-flow distal branch.
Using BRs 3 and 4, we were able to preserve the branched structure of the vascular network as well as render a form of diameter control during flow-mediated EC migration.
However, our simple bifurcation rules were based solely on randomness and arbitrarily chosen probabilities and lack a mechanism as to how and why bifurcations remain stable. Therefore, we took steps to create a bifurcation rule which contained more physiologically relevant mechanisms based on what we know about flow-mediated vascular remodelling in vivo. As mentioned previously, experiments involving noncanonical Wnt5a signalling suggest competitive interplay between flow-based polarisation/migration and collective cell junction communication. Based on these observations, we designed a mechanistic bifurcation rule (BR 5) through which the probability of each branch "attracting" incoming cells is given by a weighted average of two probability components: one due to shear stress (flow-migration term), and one due to cell number/junction force transmission (collective cell behaviour term) (see Methods, . We define the probability of an EC choosing a branch due to shear stress as that branch's contribution to the shear stress ratio (as ECs are attracted towards regions of higher shear). Likewise, probability due to junction forces is defined as the branch's contribution of the cell number ratio (as the resulting net force direction will favour the branch with larger number of cells). Use of ratios to define probability ensures that our probability definitions always rest between the required range of 0 and 1. The strength of influence each component has on EC decisions is governed by the parameter α: shear stress probability scales relative to α, and cell number probability relative to 1-α. Unlike in BR1 and BR2, ECs will choose the preferred direction of migration only probabilistically (i.e., it is still possible for an EC to choose the non-preferential direction of migration) since we want to investigate stochastic effects arising due to the relatively low number of cells in these vessel networks.
We tested this new bifurcation rule by running the model 1000 times, each with a different random seed, for different values of α (using the same 1000 random seed numbers for each value of α). Setting α to 0.0 means that only collective cell behaviour is considered at the bifurcation; 91.9% of simulations at this value of α experienced regression and bifurcation loss with no obvious preference for the proximal or distal branch (Fig 3 A, S5 File). A value of α = 1.0 means that only shear stress is considered at the bifurcation; these simulations were even more unstable, with 100% resulting in bifurcation loss shortly after 2 days of migration (Fig 3   B, S6 File). Additionally, in these simulations the low-flow distal branch was 2.7× more likely to regress than the high-flow proximal branch. Simulations at this value of α were prone to lose the bifurcation much sooner than their counterparts at α = 0.0, which didn't reach over 90% loss until 5 days of migration. It should be noted that in this bifurcation rule, branch probability continuously adapts to changes in flow and cell number (rather than being held constant as seen in BRs 3 and 4). This results in rapid adaptation of the vasculature to the various stochastic outcomes which can push the model to stable solutions that may seem unintuitive a priori. For example, in nearly 25% of cases with α = 1.0 resulted in loss of the proximal branch in favour of the longer distal path, despite the proximal path initially experiencing larger shear stress.
These results demonstrate the complex nature of the system, which can quickly and dramatically push towards different stable outcomes given accumulation of small stochastic effects. shear stress and cell number components ( for the proximal branch and for the 1 1 , n P P  2 2 , n P P  distal branch), weighted by the parameter α (Fig. 4 A). In simulations at stable values of α (e.g., = 0.45), the probability of each branch oscillated slightly over time but centred around  a relatively constant probability value (Fig. 4 B). This averaged to a higher value in the high-

Discussion
Branched vascular networks are a heavily conserved feature across animal physiology, and this geometric configuration is vital for successful transportation of blood. The fractal-like nature of embedded vascular networks provides effective transport of nutrients and waste across the tissue domain by maximising surface area available for transport, minimising diffusion distance, and reducing the energetic cost of forcing blood flow [22][23][24][25]. Therefore, during the development of these networks an optimal level of branching must be achieved.
Sprouting angiogenesis produces an initially hyperbranched network with numerous redundant flow paths, and the remodelling phase works to optimise this network by removing some branches whilst maintaining others at key locations. As of now, the exact mechanisms as to how ECs within developing networks collectively decide on what the optimal number of branches remains unknown. Franco et al. demonstrated how ECs responding to shear stress differences at vessel bifurcations is a mechanism triggering regression and loss in branching [5]. However, simulations of blood flow in developed vascular networks reveal that numerous bifurcations with shear stress differences can exist in late-stage/post-vascular remodelling [21], indicating that shear stress differences are required but not sufficient for regression and bifurcation loss. This brings us to the fundamental question of our research: what determines if a bifurcation is to be removed or preserved during flow-mediated vascular remodelling?
For a vessel to persist during remodelling, the net flux of cells entering and leaving cannot be negative. Therefore, upon reaching a vessel bifurcation, cells must finely tune migration along both possible paths if the bifurcation is to remain. Flow-migration coupling demonstrates that cells have a tendency to migrate from low to high-shear vessels [5,6,8].
However, these cells must migrate as a collective with adherens junctions intact in order to maintain fluid barrier function and prevent leakage. Therefore, the directional cues due to force transmission at adherens junctions must also play a role in influencing cell migration decisions [11,26]. However, how these migratory cues are integrated at the cellular level to promote or avoid vessel loss is unknown and currently challenging to directly prove experimentally. Our findings demonstrate that shear stress and junction force signalling interact competitively as mechanisms for determining cell migration paths at bifurcations. Furthermore, we predict that cells must integrate both cues nearly equally to maximise bifurcation stability. Dominance by one over the other leads to regression and loss of the bifurcation. The competitive nature of these two cues arises hydrodynamically from the inverse relationship between vessel diameter (and hence cell number) and resistance to flow. These findings also imply a previously unreported emergent mechanism protecting highly dynamical developmental networks against bifurcation loss: a branch that is narrowing due to a net efflux of cells will have a low junction force transmission at the bifurcation, but the decreased diameter will cause a spike in resistance and pressure drop within the vessel. This pressure spike will increase the shear stress within the vessel, attracting new cells in order to restore the vessel and prevent collapse and regression. from the distal pathways that prevents shunt formation. The authors proposed cell-cell signalling via gap junctions as a possible mechanism, but no unequivocal demonstration of this mechanism has been proposed to date. Our findings on the competitive interplay between shear sensing and junction force transmission during remodelling suggest that this additional signal that allows the distal pathways to remain intact could be force transmission through the collective endothelium during migration.
In the current study, we did not explicitly consider the molecular regulation of flowmigration coupling during angiogenic remodelling. However, several important signalling regulators of bifurcation stability have been implicated. For example, BMP-ALK1-SMAD signalling facilitates cell migration and vessel stability in low shear stress environments to prevent excessive regression [7,28]. VEGF3 modulates EC flow-migration coupling by influencing shear stress sensitivity [29]. Notch signalling facilitates polarisation against flow and artery-vein specification [30][31][32]. Lastly, signalling from vascular mural cells promotes vascular stability and survival [33], and the specific role of mural signalling during angiogenic remodelling and vessel regression is currently unknown. An advantage of our modelling approach is that we implicitly account for both flow-and collective-directed migratory cues and utilise a single parameter, α, to control the relative weight ECs place on each when selecting migration paths at bifurcations. Future work will investigate functional formulations of α that can capture this molecular regulation.
In conclusion, we have designed the first agent-based model of EC flow-migration coupling during developmental angiogenic remodelling in order to interrogate the cellular dynamics at bifurcations leading to vessel preservation or regression. We found that bifurcation stability can be achieved through a combination of flow-based and collective-based migratory behaviour, and these two factors can interact cooperatively or competitively depending on the scenario at the bifurcation. In cases of competition where EC migratory paths split, weighting both factors equally resulted in maximum stability and minimum loss in branching. Our findings were robust across numerous changes in the model, suggesting we have uncovered an inherent property of the physical system rather than an artificial construct of the model and its assumptions. Furthermore, our work provides a theoretical basis for the experimental investigation of bifurcation stability during angiogenic remodelling. Future theoretical work will investigate emerging behaviour in complex network environments (e.g., full vascular plexus models) and mathematical modelling of the molecular regulation behind flow-migration coupling. In particular, we are interested in instances in which the normal flowmediated remodelling process is disrupted, leading to pathological structural abnormalities and functional shunting of vascular beds (e.g., arterio-venous malformations).

The flow boundary-value problem
The ABM of flow-migration coupling consists of a network of blood vessels composed of migratory ECs represented by agents. This vessel network consists of a feeding vessel serving as a flow inlet which branches into a shorter proximal branch and longer distal branch.
The two branches then reconnect with each other at the draining vessel, which serves as the flow outlet. The distal branch is twice the length of the proximal branch, which is twice the length of the feeding and draining vessels (Fig 1 A). Flow, q, is driven by differences in pressure, p, along the network and pressure boundary conditions prescribed at the inlet   where w is the lateral width of an EC. Our model currently represents cells as rigid bodies with constant axial length, lateral width, and surface area. With the lumen diameter, we can then calculate the vessel conductance to flow (i.e., the amount of flow generated by a given pressure difference) as,

EC migration, boundary conditions, and bifurcation rules
The ABM then takes a migration step within the current flow environment. For the majority of segments within the network, this simply involves moving the cells it currently has to its upstream neighbour while receiving the cells from its downstream neighbour. We included a diffusion-like intercalation component in our migration step that smooths sharp fluctuations in vessel diameter. During each migration step, if the number of cells leaving a segment is greater than the number of incoming cells (and the number of incoming cells is greater than zero), then one cell was randomly chosen to stay behind and not migrate during that step. Note that this intercalation was not implemented at vessel bifurcations to avoid interference during our analysis of EC decision behaviour at these locations.
There are several exceptions to the migration step that need to be handled specifically.
The first exception is to handle cell boundary conditions, namely the case of cells migrating out at the inlet of the network as well as the condition of cells entering the network at the outlet.
We The remaining migration exceptions involve the two bifurcations within the network.
The bifurcation at the left of the network (near the inlet) is a flow-divergent bifurcation which means EC migration paths converge. At this bifurcation, incoming cells from both branches only have one choice to migrate against the flow which is to join at the feeding vessel before exciting through the inlet. The bifurcation on the right of the network (near the outlet) presents a much more interesting case. This is a flow-convergent bifurcation where EC migration paths diverge: cells approaching this bifurcation have two choices available when migrating against the flow. Due to the difference in path lengths between the proximal and distal branch, the flow and shear stress in the proximal segment (hence referred to as branch 1) is higher than the distal segment (hence referred to as branch 2).
Since the actual mechanisms regulating EC migration at vessel bifurcations during flow-mediated remodelling are unknown, we designed and implemented various "bifurcation rules" for our agents upon approaching the flow-convergent bifurcation and observed the emergent outcomes on the remodelled vessel network. In our first bifurcation rule (BR 1), shear stress is the only determinant of which branch the ECs choose, such that .
where each branch probability is constantly updated via weighted average of the probability due to shear stress and the probability due to cell number , The shear stress and cell number probabilities come from the ratio of shear stress and cell number between the two branches, respectively, , 1 2 i ni n P n n   Note that Eq 6 can be substituted into Eq 13 to express the shear stress probability as the ratio of the multiplicative combination of cell number and the pressure difference over the branch segment, .

Overview of simulations and data analysis
For a full list of all parameters and variables with the model as well as values used in the simulations, please see Table 1. We performed several simulations for each bifurcation rule and quantified the results in order to characterise the impact each rule had on network remodelling. Mean diameter of the proximal and distal branch was measured by taking the mean of the diameter across all adjacent segments in that branch. The same number of segments was included in this mean for both the proximal branch and distal branch, even though the distal branch was twice in length. The first two bifurcation rules, BR 1 and BR 2, do not depend on any randomness so only a single simulation was performed. For the remaining rules, we varied the seed number for the random number generator accordingly. We generated a list of 1000 seed numbers between 1 and 10 9 using the randint function in the random.py Python library. Each simulation utilised a different number from this list as its seed number, and the same list was used across all other parameter and rule changes. In the basic stochastic rules BR 3 and BR 4, we only ran simulations using the first 10 numbers as this was sufficient to characterise trends outside of randomness in those cases.
When assessing the parameter in the mechanistic bifurcation rule BR 5, we varied across   . In all cases, there were temporary oscillations between shear stress and cell number probability while the bifurcation remained stable, but in each case a branch probability reached maximum likelihood (i.e., equal to 1) and the bifurcation was lost.

Results
The final task of our study was to demonstrate that the stability behaviour described previously is not exclusive to our pressure-driven A branch model with periodic cell conditions. Hence, we performed similar stability analysis across different boundary conditions and vessel network geometries. The first variation of the A branch model included a prescribed flow condition at the inlet as opposed to prescribed pressure (S10 File). In this model, the pressure at the inlet varies in order to set a prescribed flow at the inlet, which was set to the same initial incoming flow in the pressure-driven model. Prescribing inlet flow instead of pressure had no effect on the stability at the flow-convergent bifurcation, as bifurcation loss vs. α over time was identical to the pressure-driven results (S11 Fig). Next, we sought to determine the impact on the EC boundary conditions on stability by holding the number of incoming cells constant, meaning the total number of ECs within the network varied in order to match this condition (as opposed to the periodic condition where the total number of ECs was held constant) (S12 File). Simulations with this condition resulted in a similar asymmetric saddle-shape while sweeping over the range of α; however, the global minimum of bifurcation loss was shallower (i.e., less stable) when compared to the periodic cell conditions (S13 Fig). In our final variation of the A branch model, we sought to determine how the initial shear stress difference at the bifurcations impacted stability. In the original A branch model, the shear stress difference at the bifurcation had an initial ratio of 2:1, proximal to distal, due to the distal path being twice the length (and hence twice the resistance) of the proximal path. We therefore created a new version of the A branch model in which the distal path was 10× the length of the proximal path, creating an initial shear stress ratio of 10:1. Stability analysis with this model resulted in a similar asymmetric saddle to the 2:1 case, although the global minimum was shallower (less stable) and shifted slightly to the left centred around α = 0.375 (S14 Fig). The A branch model is a basic representation but still includes some "network effects" as flow is split at the flow-divergent bifurcation before re-joining at the flow-convergent bifurcation. This model has some advantages, specifically in that cell behaviour at the flowconvergent bifurcation doesn't propagate to the flow inlet as cells in both paths combine at the flow-divergent bifurcation. However, varying the shear stress difference at the bifurcation in this model is not straightforward and requires a change in geometry (i.e., change in path lengths between the two branches). Thus, we created a simplified model of a flow-convergent bifurcation (hence referred to as the Y branch model) in which incoming flow from a left and right branch combine at a single bifurcation before exiting via the outlet (S15 File). Using this model, we can set the inlet pressure in both the left and right branches to be equal in order to achieve a 1:1 shear stress ratio (i.e., no shear stress difference) at the bifurcation (S16 Fig). Stability across the range of a in this model was very similar to the pressure-driven A branch model, with a global minimum across values between 0.3 and 0.6. The branch probabilities in this model initialise at 0.5 as the cell number and shear stress in both branches is the same and exhibit similar competitive oscillations which result in bifurcation stability (S17 Fig). Finally, we implemented the Y branch model with inlet flow conditions in both branches that allows us to vary the initial shear stress difference at the bifurcation at ratios of 1:1, 1:10, and 1:20 (S18 File, S19 File, S20 File). Each of these cases resulted in a similar saddle shape in the stability surface, although this saddle became shallower (less stable) and shifted to the left (towards α = 0.0) as the shear stress difference increased (S21 Fig). While there was no difference in mean diameter between the two branches with a shear ratio of 1:1, increasing the ratio between shear stress led to larger differences in diameter between the two branches (S22 Fig).

Discussion
Our stability findings were robust across changes in boundary condition, vessel geometry, and shear stress initialisation at the bifurcation, strongly suggesting that our findings are truly inherent to physical system we are representing and not a manufactured outcome of any one particular model configuration. Inlet flow conditions made no impact of stability at the flow-convergent bifurcation in the A branch model, most likely due to the face that cells from both branches recombine at the flow-divergent bifurcation rendering the inlet free from EC behaviour at the flow-convergent bifurcation. In the inlet flow version of the Y branch model, which does not include this recombining effect, we found slightly different results but the classic asymmetric saddle shape of stability remained intact. Changing the initial conditions for shear stress (and hence shear stress probability) could affect the nature of the competitive oscillations between the two components; however, when we varied these initial conditions (either by changing the resistance of the distal branch in the A branch model or the flow ratio between branches in the Y branch model) the classic asymmetric saddle shape of stability was maintained albeit more shallow and shifted to the left towards α = 0.0. With an increased initial shear stress difference, simulations that favoured junction forces resulted in fewer bifurcations lost, although results were still not as stable when compared to simulations with less of a shear stress imbalance. Changes in the boundary condition governing the number of incoming cells had the greatest impact on bifurcation stability. Simulations which held the amount of incoming cells constant exhibited a general decrease in the total number of cells within the system, and although we found the same asymmetric saddle shape the bifurcation was inherently more unstable with 40% of simulations losing the bifurcation by day 5 of migration (as compared to around 20% found with the Periodic cell condition). These findings suggest that one of the most significant regulators of bifurcation preservation may be the level of new incoming cells: a steady source of new cells greatly enhances the chance of a bifurcation remaining during remodelling. Likewise, a bifurcation which is experiencing reduced or disrupted levels of incoming cells may be more prone to instability and regression. At this point it is unclear what acts as the source of cells in the developing network in vivo. Proliferation may be isolated to the sprouting front, meaning new ECs must arrive via travelling down the network to the remodelling zone. Additionally, large vessels such as veins could also act as a reservoir of new cells. More experimental evidence is required in characterising the source of new ECs in the developing retina in order to investigate the matter further.