A Mathematical Model of Lymphangiogenesis in a Zebrafish Embryo

The lymphatic system of a vertebrate is important in health and diseases. We propose a novel mathematical model to elucidate the lymphangiogenic processes in zebrafish embryos. Specifically, we are interested in the period when lymphatic endothelial cells (LECs) exit the posterior cardinal vein and migrate to the horizontal myoseptum of a zebrafish embryo. We wonder whether vascular endothelial growth factor C (VEGFC) is a morphogen and a chemotactic factor for these LECs. The model considers the interstitial flow driving convection, the reactive transport of VEGFC, and the changing dynamics of the extracellular matrix in the embryo. Simulations of the model illustrate that VEGFC behaves very differently in diffusion and convection-dominant scenarios. In the former case, it must bind to the matrix to establish a functional morphogen gradient. In the latter case, the opposite is true and the pressure field is the key determinant of what VEGFC may do to the LECs. Degradation of collagen I, a matrix component, by matrix metallopeptidase 2 controls the spatiotemporal dynamics of VEGFC. It controls whether diffusion or convection is dominant in the embryo; it can create channels of abundant VEGFC and scarce collagen I to facilitate lymphangiogenesis; when collagen I is insufficient, VEGFC cannot influence the LECs at all. We predict that VEGFC is a morphogen for the migrating LECs, but it is not a chemotactic factor for them.


Introduction
The lymphatic system of a vertebrate plays many roles in health and in diseases. Most importantly, it drains the interstitial fluid of its tissues back to the blood vasculature, thereby maintaining tissue homoeostasis and absorbing intestinal lipids (Margaris and Black 2012;Schulte-Merker et al. 2011). Furthermore, various immune cells reside in the lymph nodes distributed throughout the lymphatic system; they filter the circulating lymph (Margaris and Black 2012). If the lymphatic system malfunctions, a medical condition called lymphoedema ensues; it is characterised by swelling and pain due to a build-up of interstitial fluid (Margaris and Black 2012).
In a vertebrate, lymphatic vessels are present in most organs except avascular tissues like cartilage (Schulte-Merker et al. 2011;Louveau et al. 2015). Like Roose and Tabor (2013), we will classify them into primary and secondary lymphatics. Margaris and Black (2012) is a detailed review of both categories. The primary lymphatics, also called initial lymphatics, are the entry points of a lymphatic system. They are lined by a monolayer of nonfenestrated lymphatic endothelial cells (LECs). They drain their surrounding tissues of excessive fluid passively, a process driven by fluctuations in their interstitial pressures. The resulting lymph is delivered into the larger secondary lymphatics. Also known as collecting ducts, they have walls that contain smooth muscle cells to propel lymphatic circulation by contractions; the muscles, arteries, and organs nearby also add to the propelling forces. The secondary lymphatics drain into various veins, thereby returning lymph to the vertebrate's blood vasculature.
How such an important and complex structure develops is incompletely understood. In this paper, we will investigate lymphangiogenesis in zebrafish (Danio rerio) embryos. Lymphangiogenesis is the development and proliferation of new lymphatics by sprouting from veins and/or any pre-existing lymphatic structures (Ji 2006). Zebrafish is a model organism widely used for studying vascular development (Gore et al. 2012). According to Florence Sabin's conceptual model (Sabin 1902), the lymphatic vasculature of a vertebrate stems from the blood vasculature. This mechanism is generally accepted by the scientific community nowadays (Schulte-Merker et al. 2011). In zebrafish, various venous origins contribute the precursor cells which will form the trunk lymphatics, the facial lymphatics, the lateral lymphatics, and the intestinal lymphatics (Koltowska et al. 2013). Our focus is on the trunk lymphatics. The developmental steps which generate the lymphatic vasculature in a zebrafish trunk are illustrated in Fig. 1 and described as follows.
Within 24 h post-fertilisation (HPF), Wnt5b secreted by the endoderm commits the cells in the ventral wall of the posterior cardinal vein (PCV) to the lymphatic fate; the resulting LECs will translocate to the dorsal side of the PCV by 30 HPF (Nicenboim et al. 2015). At 32 HPF, most of the blood vasculature is fully formed, including the dorsal aorta (DA), the PCV, a set of intersegmental arteries (aISVs), and a pair of dorsal longitudinal anastomotic vessels (DLAVs) (Koltowska et al. 2013).
At around 36 HPF, 30 pairs of secondary sprouts emerge from the dorsal side of the PCV and migrate dorsally (van Impel and Schulte-Merker 2014). The LECs constituting these sprouts only exit the PCV when they are stimulated by the growth factor VEGFC (Hogan et al. 2009). In mice at least, only the LECs that have exited the veins can express podoplanin (Koltowska et al. 2013). Although podoplanin is absent in Fig. 1 (Color figure online) Developmental steps that generate the lymphatic system in the trunk of a zebrafish embryo. a-d A slice of the trunk cut in the ventral-dorsal direction, so they depict the developmental events in the anterior-posterior view. This particular slice of the trunk has a pair of intersegmental arteries (aISVs) and a pair of lymphatic sprouts, one of which fuses with an aISV to from an intersegmental vein (vISV). There are 30 slices like this one in the trunk. When the parachordal lymphangioblasts (PLs) reach where the thoracic duct and the dorsal longitudinal lymphatic vessel lie in the ventral-dorsal slice depicted, they migrate anteriorly and posteriorly to connect with the PLs from the remaining 29 slices zebrafish (Chen et al. 2014), the genetic programmes regulating lymphangiogenesis in zebrafish and mice are more similar than different, as argued in van Impel and Schulte-Merker (2014). This leads us to assume that the PCV-derived LECs will change their gene expression profile after exiting the PCV.
At approximately 48 HPF, half of the secondary sprouts are already fused with their adjacent aISVs to form a set of intersegmental veins (vISVs); the remaining sprouts aggregate in a region named horizontal myoseptum, forming a pool of parachordal lymphangioblasts (PLs) (van Impel and Schulte-Merker 2014). The horizontal myoseptum expresses the ligand Cxcl12a which binds to the receptor Cxcr4 expressed by the LECs constituting the sprouts, thus ensuring the dorsally migrating LECs turn laterally when they reach the horizontal myoseptum . Further guidance cues for the LECs are thought to be provided by the motor neuron axons positioned along the horizontal myoseptum . After the LECs form the pool of PLs, they continue to express Cxcr4; they will migrate both ventrally and dorsally along their adjacent aISVs which express the ligand Cxcl12b .
Before 120 HPF, the PLs form the thoracic duct (TD) between the DA and the PCV, as well as the dorsal longitudinal lymphatic vessel (DLLV) below the DLAVs (van Impel and Schulte-Merker 2014). These two lymphatic vessels are connected via a set of intersegmental lymphatic vessels (ISLVs) which are close to the aISVs (van Impel and Schulte-Merker 2014). At this stage, the PCV expresses Cxcl12a and the DA expresses Cxcl12b, thus ensuring the ventrally migrating PLs will end up between the two blood vessels . Once they reach where the DLLV and TD should lie in a ventral-dorsal slice, the PLs will migrate anteriorly and posteriorly to connect with the PLs from the other ventral-dorsal slices in the trunk, thus ensuring the two lymphatic vessels are continuous.
There are several missing details in this developmental process. The LECs exit the PCV under the influence of VEGFC. During their dorsal migration, their gene expression profile probably changes, similar to their counterparts in a mouse embryo. Although we know that Cxcl12a causes the LECs to aggregate along the horizontal myoseptum, we do not know what causes them to migrate dorsally instead of ventrally or laterally from the PCV. Neither do we know what changes their gene expression profile during migration. In short, we are uncertain about what happens between (b) and (c) in Fig. 1.
We know that VEGFC promotes survival, proliferation, and migration in LECs through the PI3K/AKT and RAS/RAF/ERK signalling pathways; the PI3K/AKT pathway regulates their migration, while the RAS/RAF/ERK pathway controls lymphatic fate specification (Mäkinen et al. 2001;Deng et al. 2013). A possibility is that VEGFC is more than a growth factor for the PCV-derived LECs. It may be a chemotactic factor and a morphogen too. By chemotactic factor, we mean a chemical which directs the migrating LECs dorsally to the horizontal myoseptum. By morphogen, we mean a chemical which provides positional information to the LECs so that they alter their gene expression after exiting the PCV. In Sect. 2, we will build a mathematical model of the spatiotemporal dynamics of VEGFC in the trunk of a zebrafish embryo. In Sect. 3, we will solve the model numerically under different conditions to explore the aforementioned possibilities. In Sect. 4, we will integrate the simulation results into answers to our research questions about lymphangiogenesis.

Development of the Mathematical Model
In the following subsection, we will represent a zebrafish's trunk with a simplified geometry. Then, we will use Brinkman's equation to model the interstitial flow in the trunk. After that, we will use different forms of the reaction-diffusion-convection equation to model the reactive transport of VEGFC in the interstitial space of the trunk, as well as the changing composition of the interstitial space itself. We will complete the model by connecting the composition of the interstitial space to the interstitial flow. The resulting mathematical model will be a single framework which integrates these biochemical and biophysical phenomena. Then, we will parametrise, nondimensionalise, and simplify this mathematical model.

Geometry
According to van Impel and Schulte-Merker (2014), the blood and lymphatic vasculatures in a zebrafish trunk are spatially periodic in the anterior-posterior direction as defined in Fig. 1. Exceptions are the three sets of intersegmental vessels: the aISVs, the vISVs, and the ISLVs, which appear at certain points on the anterior-posterior axis only; they extend in the ventral-dorsal direction as defined in Fig. 1. The secondary sprouts emerge next to the aISVs, so the three sets of intersegmental vessels coalign in 30 ventral-dorsal slices of the trunk (Isogai et al. 2003). Two adjacent slices are about 75 µm apart (Coffindaffer-Wilson et al. 2011b), so they can be considered independently. Each slice is similar to the one shown in Fig. 1. We will take advantage of these features and model one slice only.
However, we will not model the aISVs because our interest is from 36 to 48 HPF, the period when the PCV-derived lymphatic progenitors migrate to the horizontal myoseptum and differentiate en route. These events are not dependent on the aISVs (Bussmann et al. 2010).
There are a pair of DLAVs, but the distance between them is small. Representing them as two separate tubes requires a high-resolution grid, so we will model one DLAV only and double the flux into this vessel.
Based on the above assumptions, we can build an idealised geometry of the trunk between 36 and 48 HPF. The geometry is shown in Fig. 2. The LEC is located halfway between the DA and the PCV. It represents an LEC on its way to the horizontal myoseptum, but it is stationary in our model. For the purpose of model development, we will divide the geometry into two domains: the LEC and the interstitial space, which is the whole geometry minus the LEC. The dimensions of the geometry and its internal structures are summarised in Table 1.

Interstitial Flow
Next, we will consider the interstitial flow in the trunk. It is driven by the pressure differences between the zebrafish's blood vasculature, interstitial space, and lymphatic vasculature (Swartz and Fleury 2007). Clearly, our representation of the zebrafish trunk does not include a lymphatic vasculature. However, the blood circulation in a   (Iida et al. 2010), so there is already an interstitial flow at the beginning of our time frame of interest. We need to incorporate this physical phenomenon into our mathematical model. On the other hand, we will not consider the flow's effects on the LEC in our model trunk. In general, the shear stresses from a flow can induce intracellular and functional changes in cells (Shi and Tarbell 2011;Ng et al. 2004). However, the intracellular details necessary for the calculation of its mechanical responses are beyond the scope of this tissue-level model. Our mathematical model of the interstitial flow relies on several assumptions. First, as in Coffindaffer-Wilson et al. (2011a), we will assume that the DA has the highest blood pressure in a zebrafish. The images in Coffindaffer-Wilson et al. (2011b) show that the DA, PCV, and DLAV have diameters comparable to a single cell. It is therefore reasonable to treat them as leaky capillaries (Jain 1987). The high DA pressure will force blood plasma into the interstitial space by paracellular transport, making the DA the inlet for fluid flow in our model. Second, we will assume a constant density for the resulting interstitial fluid. Third, we will only model the interstitial flow in the interstitial space because the LEC is separated by its cell membrane. Fourth, we will assume there are no sources or sinks of fluid in the interstitial space. Fifth, we will use a constant permeability for all three blood vessels because they are all assumed to behave like one-cell-thick capillaries. Sixth, we will ignore the pulsating nature of blood flow in this mathematical model. Finally, we will assume that the interstitial flow is at a steady state.
The interstitial space consists of the aforementioned interstitial fluid and an extracellular matrix (ECM), the latter of which is a porous medium. Therefore, the interstitial flow can be described by Darcy's law. However, Darcy's law does not permit the use of no-slip boundary conditions on the surfaces of internal structures, such as the blood vessels and the LEC in our geometry. More significantly, Darcy's law assumes a homogeneous medium. In the next subsection, we will expand the model to include the remodelling events which degrade the ECM wherein channels may form. Darcy's law cannot model these regions accurately. Brinkman's equation can overcome both limitations. Using P (mmHg) to represent the pressure field in the interstitial space, μ (cP) the dynamic viscosity of the interstitial fluid, κ (cm 2 ) the specific hydraulic conductivity of the ECM, and u (µm/s) the interstitial fluid velocity, we can write Brinkman's equation as There are two dependent variables in Eq. (1): P and u, so we need another equation to define the flow problem. Because the interstitial fluid has a constant density and there are no sources or sinks of it in the interstitial space, conservation of mass is given by To solve the flow problem, we need some boundary conditions. The fluxes out of the DA and into the PCV and the DLAV can be modelled by an equation describing the permeability of a vessel (Jain 1987). It is a linear relation between a transvascular flux and the transvascular pressure drop driving it. We will define x as the position vector in our geometry and n as a normal vector pointing out of the domain it resides in. Because the three normal vectors on the blood vessels point out of the interstitial space, they point into the vessels. Our definition also means a mass flux into the interstitial space is positive.
We will use ρ (kg/m 3 ) to represent the density of the interstitial fluid, L DA (cm/Pa/s) the DA vascular permeability, and P DA (mmHg) the pressure inside the DA. Math-ematically, the relation gives the mass flux from the DA surface (∂Ω DA ) into the interstitial space as We can derive the boundary conditions on the PCV and DLAV surfaces (∂Ω PCV and ∂Ω DLAV ) along the same line to give The multiplicative factor of 2 in Eq. (5) is there because we are representing two paired DLAVs as one vessel.
Finally, we will impose no-slip boundary conditions on the four outer boundaries of the geometry, collectively labelled ∂Ω x,y , and the LEC surface, seen from the interstitial space domain, ∂Ω LEC/IS+ . They are represented by

Reactive Transport of VEGFC and Extracellular Matrix Remodelling
In this subsection, we will add a biochemical reaction network to the mathematical model. We will model the transport phenomena of the participating biochemical species too. VEGFC is synthesised as a preproprotein called proVEGFC; it has an N-terminal signal sequence followed by an N-terminal propeptide, then the VEGF homology, and finally a cysteine-rich C-terminal segment (Joukov et al. 1996(Joukov et al. , 1997Siegfried et al. 2003). proVEGFC undergoes cleavage intracellularly and extracellularly (Joukov et al. 1997). After intracellular processing, proVEGFC will become a tetramer which has a molecular weight of 120 kDa and it will be secreted (Joukov et al. 1997). The secreted tetramer will bind to a VEGFR3 receptor on an LEC. On the cell surface, it is cooperatively cleaved by CCBE1 and ADAMTS3 (Jeltsch et al. 2014). Our investigation is concerned with the spatiotemporal dynamics of VEGFC on the tissue level, so we are not interested in these events which occur on the cellular level. Therefore, we will not model any cleavage events of VEGFC, intracellular or extracellular, and VEGFC-VEGFR3 binding. It follows that VEGFC denotes the tetramer only in this investigation and it is limited to the interstitial space domain.
We have not discussed the properties of the ECM yet. Its major structural components include different kinds of collagens and glycosaminoglycans (Lutter and Makinen 2014). Collagens make up more than two-thirds of the ECM protein content in many soft tissues (Swartz and Fleury 2007). According to Prockop and Kivirikko (1995), collagen type I is the most abundant protein in humans. More specifically for our study, LECs are mainly surrounded by fibrillar type I collagen in general (Wiig et al. 2010;Paupert et al. 2011). The ECM of an embryo regulates its lymphangiogenic processes in several way (Lutter and Makinen 2014). First, it confers structural support and stability to the embedded cells, tissues, and organs, but it is also a barrier to cell migration. Second, the ECM contains components that can bind to a myriad of cell surface receptors, thus inducing intracellular changes. Third, the ECM can bind to growth factors, thus sequestering them and creating concentration gradients. It is the third function that interests us in this investigation. We will model the transport of VEGFC in the interstitial space domain where it interacts with the ECM. Since LECs are generally surrounded by fibrillar type I collagen, we will treat the ECM as pure collagen I in our model. VEGFC binds to heparan sulphate (Lutter and Makinen 2014), but we do not know whether it binds to collagen I. In order to mimic VEGFC's interactions with the ECM without modelling heparan sulphate explicitly, we will assume that VEGFC binds to collagen I reversibly in an 1:1 stoichiometric ratio.
An ECM is not inert and undergoes constant remodelling. According to Helm et al. (2007), LECs secrete a protease called matrix metallopeptidase 9 (MMP9) to degrade collagen, thereby rendering their surrounding ECM more conducive to their migration. According to Bruyère et al. (2008), LECs can produce and activate another protease called matrix metallopeptidase 2 (MMP2) to regulate lymphangiogenesis. Commenting on Bruyère et al. (2008), it is argued in Detry et al. (2012) that MMP2 is more important than MMP9. This theory explains lymphangiogenesis in terms of LEC migration through an interstitial collagen I barrier and a collagenolytic pathway driven by MMP2 (Detry et al. 2012). In this investigation, we will consider the production and activation of MMP2 in the LEC domain, as well as the degradation of collagen I by MMP2 in the interstitial space domain. A conceptual model of these MMP2related events is proposed in Karagiannis and Popel (2004). In this conceptual model, proMMP2, TIMP2, and MT1-MMP act cooperatively to activate proMMP2 to form MMP2. Although MT1-MMP is restricted to the surfaces of LECs, we will disperse the MT1-MMP molecules uniformly in our LEC domain to simplify the mathematics. In our model, the cooperative action occurs in the LEC domain to produce the mature MMP2. However, proMMP2, MMP2, and TIMP2 can all diffuse into the interstitial space domain. In the interstitial space domain, TIMP2 can bind to and inhibit MMP2 reversibly. Karagiannis and Popel (2006) is a mathematical modelling study based on this conceptual model and is an inspiration for our study. It is possible for the interstitial flow to affect the ECM's composition, either mechanically or by stimulating the LEC to produce or degrade ECM components. Nonetheless, we will assume that the ECM's behaviour is dominated by collagen I and MMP2 dynamics.
Combining the biochemical events described in this subsection, we can construct the overall biochemical reaction network shown in Fig. 3. We will assume that the number of MT1-MMP molecules in the LEC domain is constant, meaning its production rate equals its shedding rate. This assumption allows us to ignore shedding in this study because the shedded species do not interact with the modelled species. MT1-MMP has its own collagenolytic activity too, but it is localised to the LEC domain. In our model, the LEC is a stationary circle devoid of collagen I, so we do not need to model the collagenolytic action of MT1-MMP. Finally, we will not model degraded collagen I explicitly. As far as we are aware, degraded collagen I does not affect lymphangiogenesis. We will use a set of reaction-diffusion-convection equations to model the spatiotemporal dynamics of the mobile species in the interstitial space. We will use C i (M) to represent the molar concentration of species i; t (s), time; D eff i (µm 2 /s), the effective diffusivity of species i; ω, the volume fraction where diffusion occurs; u (µm/s), the velocity from our mathematical model of the interstitial flow; R IS i (M/s), the net rate of production of species i at a point in the interstitial space. The equation Diffusion does not occur in collagen I fibrils or the fluid associated with them. The volume into which a molecular species can diffuse should be based on the specific 'wet' weight of collagen I. Denoting the partial specific volume of hydrated collagen I by v C1h (cm 3 /g) and the combined mass concentration of free and VEGFC-bound collagen I by [Cl] m (kg/dm 3 ), we can use a relation from Levick (1987) for ω, The effective diffusivity of each diffusible species can be calculated by Ogston's equation (Ogston et al. 1973). Labelling the diffusivity of species i in pure interstitial fluid by D ∞ i (µm 2 /s), the volume fraction of dry collagen I fibrils (without associated water molecules) by ν, the radius of a collagen I fibril by r f (µm), and the Stokes-Einstein radius of species i by r s,i (µm), the equation is Using v C1 (cm 3 /g) to label the partial specific volume of dry collagen I, the equation for the dry volume fraction of collagen I is given in Levick (1987), Denoting the Boltzmann constant by k B (1.380648813 × 10 −23 J K −1 ) and temperature by T (K), the Stokes-Einstein radius is given in Einstein (1905) as We will use a set of ordinary differential equations to model the temporal dynamics of the immobile species in the interstitial space, In the LEC domain, there is neither collagen I nor an interstitial flow. With R LEC i (M/s) being the net rate of production of species i at a point in the LEC, the governing equations for the mobile species there are reaction-diffusion equations of the form In the LEC domain again, the governing equations for the immobile species are a set of ordinary differential equations as follows, The reaction terms, R IS i and R LEC i , can be worked out from the biochemical reaction network in Fig. 3. We will apply mass action kinetics to most reactions. Mass action kinetics assumes that the rate of a reaction is the product of a rate constant and the concentrations of the participating reactants. An exception is the enzymatic degradation of collagen I by MMP2, to which we will employ Michaelis-Menten kinetics. If we assume the complex, MMP2·collagen I, is at a steady state, we can approximate collagenolysis by MMP2 as a single second-order kinetic process. This approximation is used and justified in Karagiannis and Popel (2004). We can do the same for the last two steps in the activation of MMP2 by assuming a steady state for the quaternary complex. This approximation is also used and justified in Karagiannis M2 and T2 complex reversibly. VEGFC binds to C1 reversibly. M2 degrades C1 catalytically. M2P, M2, T2, M2 · T2, and VEGFC degrade in the interstitial space.  (2004). There are also some 'auxiliary' reactions. proMMP2 and TIMP2 are produced in the LEC. The molecules in the interstitial space are subjected to attacks from enzymes (Gutfreund 1993). Exceptions are collagen I and VEGFC sequestered by collagen I. We are already modelling collagen I degradation by MMP2; sequestration by collagen I protects VEGFC from enzymatic attacks. The reaction terms in the interstitial space, R IS i , are given in Table 2. The reaction terms in the LEC, R LEC i , are given in Table 3.
On the four outer boundaries of the geometry and the vessel surfaces, we will impose no-flux boundary conditions for each mobile species. Our choices are justified as follows. First, there is no evidence that zebrafish lose the modelled molecules through their skin. While we could model this hypothetical mechanism, our model will not be less descriptive without it. Second, to the best of our knowledge, there is no evidence that the modelled molecules can enter the blood vessels. If they do pass through the vessel surfaces, we have no idea how much is filtered by the lining cells which may have binding receptors, for example. Anyway, they will simply degrade inside the vessels, so we will not model this hypothetical and poorly defined phenomenon. An exception is the flux of VEGFC from the DA surface into the interstitial space. VEGFC is expressed in the hypochord, the dorsal aorta, and the ventral mesenchyme of a zebrafish at 48 HPF (Hogan et al. 2009). The high pressure in the DA means convection is most significant around it. We will place the source of VEGFC on the DA's surface. This arrangement will retain the essences of VEGFC transport in zebrafish embryos and keep the model simple simultaneously. We should remind ourselves that VEGFC is not released from the blood inside the DA. Its production by a part of the DA's wall is independent of that elsewhere on the DA, so its release rate does not change along the DA.
On the PCV and the DLAV, they are On the DA, the boundary conditions for the mobile species except VEGFC are For VEGFC, with R DA VC (mol µm −2 s −1 ) being the release rate of VEGFC on the surface of the DA, the constant flux is On the LEC's surface, we will impose continuity conditions on the species which can cross the boundary between the LEC and interstitial space domains. We will denote the LEC surface seen from the interstitial space by ∂Ω LEC/IS+ . Seen from the LEC domain, the surface is ∂Ω LEC/IS− . Continuity applies to proMMP2, MMP2, and TIMP2 as follows, (20) On the LEC surface seen from the interstitial space, we will apply no-flux conditions on MMP2·TIMP2 and VEGFC, We also need a set of initial concentrations. In the LEC domain, we will assume that only MT1-MMP is present at t = 0; in the interstitial space, we will assume that only collagen I is present at t = 0. We will label the initial concentrations of these species by C MT1,0 and C C1,0 , respectively.

Connection of Extracellular Matrix Remodelling to Interstitial Flow
The interstitial flow is dependent on the composition of the interstitial space. Specifically, it depends on the dynamics of the ECM. This is obvious from Eq. (1), where κ is the specific hydraulic conductivity of the ECM. As the ECM remodels, this conductivity will also change, thereby affecting the interstitial flow.
We will define κ (cm 4 s −1 dyn −1 ) as the hydraulic conductivity of the ECM and [Collagen I] as the mass fraction of free and VEGFC-bound collagen I in the interstitial space. Using the experimental data presented in Levick (1987), we can relate the two by log κ = −2.70 log Collagen I − 14.18.
With M C1 (kg mol −1 ) being the molar mass of collagen I, [Collagen I] is given by Equation (24) assumes a density of 1 kg dm −3 for the interstitial space. It is also assumed that the combined mass of the interstitial fluid and the ECM in the interstitial space is conserved. When collagen I degrades, the products remain in the interstitial space, so the mass of a region is fixed at 1 kg dm −3 multiplied by the volume of the region. The experiments cited in Levick (1987) were carried out using a reference fluid with a dynamic viscosity of 1 cP. We can therefore convert the hydraulic conductivity to the specific hydraulic conductivity by κ = κ × 10 −2 dyn s cm −2 .

Parametrisation
We will begin our parametrisation with the interstitial flow component of the model. We need the pressures inside the three blood vessels, their permeabilities, as well as the properties of the interstitial fluid and the ECM. Table 4 summarises these parameters whose origins are explained below. In Hu et al. (2000), the DA peak systolic and end-diastolic pressures of zebrafish embryos are related to their wet body weights. At 48 HPF, the wet body weight is 0.714 mg. Ignoring the pulsating nature of blood flow, we will average the measured peak systolic (0.2433 mmHg) and end-diastolic (0.1255 mmHg) pressures corresponding to this wet body weight. This gives an estimated DA pressure of 0.1844 mmHg relative to the pressure outside the embryo. The PCV and the DLAV are parts of a closed and pumped blood circulatory system, so they must have higher pressures than the embryo's unpumped surrounding. However, in the absence of data about pressure drops in the circulatory system, back-of-the-envelope estimates are unlikely to be accurate. For the sake of simplicity, we will just set the PCV and DLAV pressures to zero. The DA pressure must drop along the vessel from the zebrafish's heart. Our geometry includes a slice of the DA only, so the pressure drop does not show up in the model. It is not important either because the lymphangiogenic processes in Fig. 1 do not occur along the anterior-posterior axis. Diffusivity of M2 0.85 × 10 −6 cm 2 /s Karagiannis and Popel (2006) D ∞

T2
Diffusivity of T2 1.10 × 10 −6 cm 2 /s Karagiannis and Popel (2006) r f Radius of a C1 fibril 2 nm Karagiannis and Popel (2006) v C1 Specific volume of dry C1 0.75 cm 3 /s Levick (1987)  We have already decided to use one vascular permeability for the three blood vessels. We will rely on Jain (1987) for this parameter. The reported measurements are concerned with species ranging from Guinea pigs to frogs. Although zebrafish is not among these species, we can use the permeability for frog skeletal muscles, 7.2 × 10 6 cm s −1 cmH 2 O −1 . The justification is that both frogs and zebrafish are cold-blooded.
The specific hydraulic conductivity of the ECM is given by the relation in Sect. 2.4. However, we need M C1 in Eq. (24). The model in Karagiannis and Popel (2006) uses a molecular weight of around 300 kDa for a collagen I fibril, which is equivalent to a molar mass of 300 kg mol −1 .
The interstitial fluid contains roughly 40% of the protein concentration of blood plasma (Swartz and Fleury 2007). Their similarities in composition allow us to use the parameter values for blood plasma. At 37 • C, the dynamic viscosity of blood plasma is 1.2 cP (Swartz and Fleury 2007) and its density is 1025 kg/m 3 (Frcitas 1998).
We will turn our attention to the reaction-diffusion-convection equation and its reduced forms next. Their parameters can be categorised into transport and kinetic parameters.
In order to calculate the volume fraction where diffusion occurs using Eq. (8) and the effective diffusivity of a species using Eq. (9), we need several parameters. We will assume a temperature of 298 K. A collagen I fibril is approximately 300 nm long and 4 nm in diameter, so r f is 2 nm (Karagiannis and Popel 2006). We can also find the required D ∞ i values in Karagiannis and Popel (2006) and Berk et al. (1993). The partial specific volume of dry collagen I and that of hydrated collagen I, v C1 and v C1h , are 0.75 and 1.89 cm 3 /g (Levick 1987). These parameters are summarised in Table 5.
The majority of our kinetic parameters about ECM remodelling are from Karagiannis and Popel (2006,2004). Their sources are various experimental studies. We are unaware of any data on VEGFC-collagen I interactions. In fact, we do not even know whether VEGFC binds to collagen I. We must rely on other data. In Köhn-Luque et al. (2013), there are reports of experimentally estimated parameters on the interactions  between VEGF (related to but different from VEGFC) and various ECM molecules like fibronectin and heparan sulphate proteoglycans. We will use the general degradation rate constant used in Hashambhoy et al. (2011). For the production rates of proMMP2 and TIMP2, we will use the secretion rates of general MMP and TIMP by endothelial cells, used in Vempati et al. (2010). These estimates are in molecules/cell/h, but we can convert them to M s −1 using the known LEC diameter of 10 µm. We are unaware of any data on the production rate of VEGFC. Therefore, we will use the secretion rate of VEGF by endothelial cells, estimated in Hashambhoy et al. (2011). These parameters are summarised in Table 6.
Finally, we need the initial concentration of MT1-MMP in the LEC and that of collagen I in the interstitial space. C MT1,0 is unavailable for endothelial cells. However, a value based on other cell types, 180,000 molecules/cell or 5.71 × 10 −7 M, is in Karagiannis and Popel (2006). In adult tissues, the concentration of collagen I ranges from 1.76 × 10 −4 to 5.29 × 10 −4 M (Levick 1987;Karagiannis and Popel 2006). We will use the midpoint of this range, so C C1,0 is 3.50 × 10 −4 M for adult tissues. The collagen content of frog embryos and larvae ranges from 4.51×10 −7 M to 2.73×10 −6 M (Edds Jr 1958). We will use the midpoint of this range, 1.59 × 10 −6 M.

Nondimensionalisation
Nondimensionalisation reduces the number of parameters in a mathematical model, gives us insights into the model in terms of its key parameters and characteristic properties, and identifies any mathematical techniques for approximation like limiting cases. We need the characteristic scales of the model in order to nondimensionalise it. To obtain the nondimensionalised spatial coordinates,x = x L , we need the length scale, L; we will use the largest dimension of the geometry, 434 µm. We will use P DA = 0.1844 mmHg to scale the pressure field,P = P P DA . Because we are interested in the period from 36 to 48 HPF, we will use a time scale, τ , of 12h. The nondimensionalised time is thereforet = t τ . Velocity is scaled likeũ = u U . The concentrations are nondimensionalised likeC i = C i C i,s . Since we are not modelling collagen I synthesis, its initial concentration is also its highest possible concentration. Therefore, we will use the initial concentration of collagen I as its scale, C C1,s = C C1,0 . We will use the adult value, C C1,0 = 3.50 × 10 −4 M. The concentrations of MT1-MMP and its two complexes will always add up to the initial concentration of MT1-MMP, C MT1,0 = 5.71 × 10 −7 M. We will use this value for C MT1,s , C MT1·T2,s , and C MT1·T2·M2P,s . We will determine the velocity scale, U (µm/s), and the remaining concentration scales while nondimensionalising the model. The scales are summarised in Table 7. Nondimensionalising the model by these scales, we will obtain a model parametrised by the dimensionless groups in Tables 8 and 9. Below are the details of how nondimensionalisation is carried out for our mathematical model.
Second, we need to nondimensionalise the equations governing the concentration fields in the interstitial space, (7) and (12). To do so, we need to introduceD eff Then, the nondimensionalised reaction-diffusionconvection equation takes the form We will define the dimensionless parameters λ 5 = v C1h M C1 C C1,s and λ 6 = v C1h M C1 C VC·C1,s . Then, we can write ω as D eff i can be expressed in terms of the dimensionless parameters λ 1,i =

s . This expression is
For the immobile species in the interstitial space, the reaction-diffusion-convection equation is reduced to Third, we will nondimensionalise the equations governing the concentration fields in the LEC domain, (13) and (14). After introducingR LEC For the immobile species in the LEC domain, the reaction-diffusion equation is reduced to Fourth, we will nondimensionalise the boundary and initial conditions of the concentrations. Equations (15-17), which apply to proMMP2, MMP2, TIMP2, MMP2·TIMP2, and VEGFC, have the following nondimensionalised forms, The boundary conditions for proMMP2, MMP2, TIMP2, and MMP2·TIMP2 on the DA's surface, represented by Eq. (18), have the following nondimensionalised form, With λ DA VC = R DA VC τ C VC,s L being the dimensionless production rate of VEGFC, we can scale C VC by setting λ DA VC = 1, so C VC,s = 1.64×10 −10 M. Then, the influx of VEGFC from the DA surface into the interstitial space, Eq. (19), has the nondimensionalised form The boundary conditions on the LEC's surface, Eqs. (20-22), are then nondimensionalised to give The continuity condition defined by Eqs. (43) and (44) applies to proMMP2, MMP2, and TIMP2. The no-flux condition defined by Eq. (45) applies to MMP2·TIMP2 and VEGFC.
The initial concentrations of MT1-MMP in the LEC domain and collagen I in the interstitial space domain areC MT1 = 1 andC C1 = 1, respectively.
Finally, we need to define the nondimensionalised reaction terms, s , in terms of nondimensionalised concentrations. We need the remaining concentration scales. We will scale C M2P by setting the dimensionless production term of proMMP2 to unity, P M2P τ C M2P,s = 1. This leads to C M2P,s = 1.14×10 −3 M. We will set the dimensionless production term of MMP2 to unity, This gives the concentration scale for MMP2, C M2,s = 3.94 × 10 −5 M. Setting the dimensionless production term of TIMP2 to unity, P T2 τ C T2,s = 1, we can obtain the concentration scale C T2,s = 6.65 × 10 −6 M. If the binding term of MMP2 and TIMP2 is set to unity,  Table 10.

Simplification
After nondimensionalisation, we can spot several opportunities for simplifying the model. First, the nonlinear term in Eq. (26),ũ (C C1 +η 2CVC·C1 ) α , can be linearised ifC C1 and C VC·C1 differ by orders of magnitude. After factorisation, the nonlinear term becomes ) −α , which can be expanded as a Taylor series at η 2CVC·C1 This expansion assumes that η 2CVC·C1 C C1 1 is true in the interstitial space domain.
The binding rate constant of VEGFC and collagen I inR IS VC·C1 is λ 31 = 1; the unbinding rate constant of the same process is λ 32 = 1.56 × 10 2 . The initial concentration of collagen I is unity, but it is zero for VEGFC bound by collagen I. It follows that C VC·C1 must be orders of magnitude smaller thanC C1 all the time. We also know that η 2 = 0.255. Taken together, the required assumption is a reasonable one to make. By the same argument, the terms beyond unity in the square brackets can be neglected, leading to a simplified version of Eq. (26), To make it easier to solve the model numerically, we will divide Eq. (47) by η 3 and define the simulated pressure,P s =P η 3 , resulting iñ Second, the blood vessels are very leaky. In the equations representing the transvascular fluxes of interstitial fluid, (28-30), the dimensionless vascular permeabilities, η DA = 1.317 × 10 14 , η PCV = 1.317 × 10 14 , and η DLAV = 2.634 × 10 14 , are all orders of magnitude larger than unity. In physiological terms, the blood vessels are so leaky that we can ignore any transvascular pressure drops. As a result, the fluxes can be replaced by constant pressures, leading tõ Third, we can linearise the diffusive term in Eq. (32),D eff i∇ (C i ω ). We will start with . We can apply (1+x) −1 ≈ 1−x +x 2 −x 3 +x 4 −· · · around x = 0 (Abramowitz and Stegun 1964) to obtain This approximation assumes that λ 5CC1 + λ 6CVC·C1 1 holds in the interstitial space. Since λ 5 = 1.98 × 10 −1 andC C1 ≤ 1 are true, λ 5CC1 1 must hold too. With respect to the second term, we know that λ 6 = 5.06 × 10 −2 andC VC·C1 C C1 , so it is orders of magnitude less than unity. The assumption is a valid one. Also, the terms beyond λ 5CC1 are negligible. The result is 1 ω ≈ 1 + λ 5CC1 . The effective diffusivity can be simplified along similar lines. The Taylor series expansion of a general exponential function around x = 0 is e x ≈ 1 + x + x 2 2! + x 3 3! + · · · We can apply it to the exponential term inD eff As usual, we need the magnitude of the exponent to be much smaller than one for Taylor expansion to work. Since λ 3 = 7.88 × 10 −2 andC C1 ≤ 1 are true, the first term in the square root is much smaller than one; λ 4 = 2.01×10 −2 andC VC·C1 C C1 mean the same for the second term. The square root λ 3CC1 + λ 4CVC·C1 is on the order of 1 × 10 −1 or smaller. Since λ 2,i is on the order of unity, this expansion is justified. For the same reason, the terms beyond −λ 2,i λ 3CC1 + λ 4CVC·C1 in the expansion are negligible, thus reducing the diffusivity to λ 1,i 1 − λ 2,i λ 3CC1 + λ 4CVC·C1 .

Numerical Experiments
In this section, we will solve the mathematical model numerically, including numerical experiments which solve the model with modified model structures and parameters. We will start by explaining how we will solve the model using COMSOL Multiphysics, including convergence studies to determine the appropriate mesh and time step sizes. Then, we will investigate the spatiotemporal dynamics of VEGFC, MMP2, and collagen I under different circumstances.

COMSOL Multiphysics Settings
COMSOL Multiphysics version 5.2 is a software package which uses the finite element method to solve partial differential equations numerically. First, the interstitial flow equations, (27) and (48), will be solved alone for the initial velocity field. Second, we will solve the time-dependent equations (32) and (35-37) together with Eqs. (27) and (48).
We will adopt a fully coupled approach to solve the model. It means that Eqs. (27) and (48) will be solved together in the first step, while Eqs. (27), (48), (32), and (35-37) will be solved together in the second step. Since our equations are all nonlinear, discretisation will result in a system of nonlinear algebraic equations. In the first step, we will use an 'automatic highly nonlinear (Newton)' solver; the second step will be tackled by the 'constant (Newton)' solver. The former has a minimum damping factor of 1 × 10 −8 ; its damping factor cannot change more than tenfold in one iteration; it terminates when the estimated relative error is less than 0.001. The latter has a constant damping factor of 0.9 and terminates when the estimated relative error is less than 0.01.
The second step runs fromt = 0 tot = 1. We will use the BDF (backward differentiation formula) method to determine the time steps adaptively. At each time point, the solutions at the previous one or two time steps are used to estimate the time derivatives at the next time step, while the stability of the derivatives determines the step size. As a starting point, we will set a maximum step size of 0.02. We will optimise this criterion after optimising our mesh setting. The mobile species are absent from the geometry initially. They need time to permeate it, so their concentrations will change drastically in the first few time steps. These time steps require tight control. Based on U = 1.371 × 10 −4 µm/s, it takes a particle 7.29 × 10 3 s to transverse 1 µm by convection. From D ∞ T2 = 1.10 × 10 −6 cm 2 /s, the highest diffusivity in our model, we know that the equivalent time is 9.09 × 10 −3 s for diffusion. Nondimensionalising the shorter travelling time, 9.09 × 10 −3 s becomes 2.10 × 10 −7 . We will therefore set the initial time step at 1 × 10 −7 .
We will perform a convergence study to determine the appropriate mesh size. In COMSOL Multiphysics version 5.2, there are several predefined mesh settings. We will solve the mathematical model with the 'fine', 'finer', 'extra fine', and 'extremely fine' mesh settings. We will run the simulations on a desktop computer with an Intel(R) Core(TM) i5-3570 CPU at 3.40 GHz and 16 GB of RAM. Figure 4 is the resulting convergence plot. The interstitial flow ranges from 0.1 to 2 µm/s in speed (Swartz and Fleury 2007). The second data point in Fig. 4 is 0.0173 µm/s, so convergence is achieved after one refinement. Thus, the 'finer' mesh setting is sufficient to capture the relevant biological information of the model. With this setting, there are 33562 elements ranging from 1.98 × 10 −5 to 0.037 in nondimensionalised length.
We will next optimise the maximum size of a time step. The equations governing u, (27) and (48), are static. A variable controlled by Eq. (35), an ordinary differential equation in time, is more suitable for this convergence study. When the maximum time step is 0.02 and the 'finer' mesh setting is used, the minimumC C1 att = 1 is 0.9991. When the maximum time step is halved to 0.01, the minimumC C1 att = 1 is 0.99909. The difference is 0.001%. The biologically relevant range ofC C1 is from 0.01 to 1 (Edds Jr 1958;Levick 1987;Karagiannis and Popel 2006). We can safely assume that a maximum time step of 0.02 is appropriate.
Our geometry is symmetric in the lateral direction. The line of symmetry runs vertically through the whole geometry and passes the centres of the blood vessels and the LEC. We will simulate half of the geometry only to save computational cost. To do so in COMSOL Multiphysics version 5.2, we will remove everything to the right of the line of symmetry from the geometry. Then, we will set the normal component of the velocity and those of the molecular fluxes to zero at the symmetry boundary.
All the simulations discussed in this paper were carried out on the aforementioned desktop computer. With the settings described in this subsection, it took 1300 s to solve the mathematical model. Henceforth, we will call this solution the primary simulation. The other simulations involved similar computational costs.

Diffusion and Sequestration Act Together
This subsection considers the primary simulation. First, we will determine the dominant transport phenomenon. Péclet number, Pe, measures the relative importance of convection and diffusion in a transport process. For our mathematical model, it is U L λ 1,i . This expression is based on the characteristic velocity scale, but the velocity magnitude varies in space. A better measure is |ũ|Pe, the maximum of which is 0.14909 att = 1 for VEGFC. Of all the diffusible species in our model, VEGFC has the smallest λ 1,i , so |ũ|Pe can only be smaller for another species. Diffusion is the dominant mode of transport in the primary simulation.
The spatiotemporal dynamics ofC VC are shown in Fig. 5. Throughout our time frame of interest,C VC peaks at the DA and decreases away from it. This symmetric distribution around the source is consistent with diffusion being the dominant mode of transport. As time passes, VEGFC achieves a wider reach by diffusion. It means that the transport of VEGFC does not equilibrate on this time scale. Because VEGFC is constantly produced, the baseline concentration of it increases with time.
Based on the simulation results, VEGFC is unlikely to be a chemotactic factor because it increases from the PCV to the DA and then decreases to the horizontal myoseptum. If VEGFC is a chemoattractant, it can guide the LEC to the DA, but the cell will remain there instead of moving further to the horizontal myoseptum. If it is a chemorepellent, the LEC should not migrate dorsally at all.
On the other hand, we argue that VEGFC is a morphogen for the LEC. The LEC is in a VEGFC gradient which increases towards the DA for the entirety of the simulation. Att = 1,C VC rises from around 0.0001 at the dorsal end of the PCV to 0.0034 at the ventral end of the DA, roughly a thirtyfold increase over 25 µm. It is estimated in Gurdon and Bourillot (2001) that a cell can read a threefold change in a morphogen's concentration over 30 µm. This estimate is based on the responses of genes to Dpp and activin gradients. The concentration profile of VEGFC in the primary simulation therefore allows it to be a morphogen. On the other hand, the viable range of known morphogens is from 1 × 10 −9 to 1 × 10 −11 M (Gurdon and Bourillot 2001).C VC in the primary simulation is on the order of 1 × 10 −13 M. Assuming sequestered VEGFC cannot function like free VEGFC, there is insufficient VEGFC. However, our estimated VEGFC production rate and VEGFC-collagen I binding rate constants are crude. If we increase λ DA VC tenfold in a numerical experiment, the baseline ofC VC will increase tenfold, but its profile's shape will remain unchanged. This is illustrated in Fig. 6a.
We will perform another numerical experiment by setting λ 25 = 0, λ 26 = 0, λ 29 = 0, λ 30 = 0, λ 31 = 0, and λ 32 = 0. Biophysically, this means we will turn off any interactions between VEGFC and collagen I. The simulated dynamics in this numerical experiment are shown in Fig. 6b. Although the peak ofC VC remains at the DA and its baseline still increases with time, the gradients are now too flat for VEGFC to be a morphogen. Sequestration by the ECM in a zebrafish embryo is necessary for VEGFC to function as a morphogen when the embryo is diffusion-dominant.
To conclude, in a diffusion-dominant zebrafish embryo, VEGFC may act as a morphogen, but it is not a chemotactic factor. For it to be a morphogen, it must bind to the ECM to steepen its concentration gradients. Sufficient VEGFC must be produced too because its production rate controls its concentration baseline.

MMP2 Acts Globally
We will now return to the primary simulation. Figure 7 summarises the simulated behaviour of MMP2. As shown in Fig. 7a, MT1-MMP gets depleted very quickly in the LEC. This rapid depletion of MT1-MMP means there is an initial burst of MMP2 production followed by a prolonged period of low production; MMP2 activation requires MT1-MMP to break up MT1-MMP·TIMP2·proMMP2. The dynamics of MMP2 production are illustrated by Fig. 7b. From Fig. 7c, we can see that MMP2 is almost homogeneously distributed. Therefore, we can consider its diffusion to be in equilibrium on this time scale. Since MMP2 degrades collagen I, the latter has a nearly homogeneous distribution too. This can be seen in Fig. 7d. The extent of this degradation is less than 0.1% and therefore negligible.
However, we know that the production rates of MT1-MMP and TIMP2 are dependent on the expressing cells' environment (Noel and Sounni 2013;Sternlicht and Werb 2001). An LEC may be able to produce more MT1-MMP or less TIMP2 when the existing stock of MT1-MMP in its environment is depleted. Furthermore, we are only modelling one LEC when each secondary sprout consists of multiple LECs. Our model certainly underestimates the collagenolytic effects of MMP2. We will perform a numerical experiment which replaces the activation mechanism of MMP2 with a constant production rate. This simplified model ignores proMMP2, TIMP2, MT1-MMP, and their complexes. Only the equations governing the interstitial flow,C M2 ,C VC , C C1 , andC VC·C1 remain. Furthermore, we will useR LEC M2 = 10 andR IS M2 = −λ 15CM2 in the modified model. In Fig. 7b, the maximumR LEC M2 is around 0.2, soR LEC M2 = 10 is a significant increase. The results of this numerical experiment are presented in Fig. 8. Similar to the primary simulation, MMP2 is almost homogeneously distributed. The flat MMP2 gradients explain why collagen I degradation does not vary significantly in space. On the other hand, the extent of degradation is significant, almost 20%. To conclude, MMP2 diffusion in a zebrafish trunk can be considered to be in equilibrium on the time scale of lymphangiogenesis. When diffusion is the dominant mode of transport, MMP2 will permeate the entire embryo and degrade collagen I almost homogeneously. This is in agreement with the conceptual model proposed by Karagiannis and Popel (2006). The diffusible protease MMP2 remodels the embryo's ECM, making it more conducive to cell migration; the LEC-bound protease MT1- Another conclusion is about TIMP2. In the primary simulation, MMP2-TIMP2 binding does not affect the homogeneous distribution of MMP2, so TIMP2 diffusion is also in equilibrium . It only changes the baseline concentration of MMP2, not its spatial distribution. In this paper, MMP2-TIMP2 binding is only considered in the primary simulation.

Convection and Asymmetry
The biologically relevant range of C C1 is from 1.59 × 10 −6 M to 3.50 × 10 −4 M. Our choice of C C1,s is 3.50 × 10 −4 M, so the initial conditionC C1 = 1 falls at the upper end of the range.
We will carry out a numerical experiment where the initialC C1 is uniformly reduced to 0.1. We will ignore collagen I degradation by MMP2 and MMP2-TIMP2 binding in the modified model, meaning it will only retain the interstitial flow Eqs. (27) and (48), as well as the governing equations for free VEGFC, collagen I, and VEGFC bound to collagen I. The two neglected phenomena can only reduce the baseline ofC C1 in a diffusion-dominant zebrafish embryo. Our manual reduction in the ini-tialC C1 has the same effect. A reduction inC C1 will affect both the transport and kinetic terms in our model. On the other hand, there are ECM components other than collagen I in a real zebrafish embryo; changing the amount of collagen I present will affect the transport properties, but the ability of components like heparan sulphate to sequester VEGFC will not be affected. For this reason, we will modify the model to control the kinetic properties. Since we are loweringC C1 uniformly by an order of magnitude, we will increase the binding terms of VEGFC and collagen I by an order of magnitude. It means the reaction terms will becomeR IS Then, we will repeat the numerical experiment after switching off the interactions between VEGFC and collagen I. This means we will further simplify the model by eliminating the governing equations for collagen I and bound VEGFC. This modified model will only have the interstitial flow Eqs. (27) and (48) in addition to the reactiondiffusion-convection equation governingC VC . The reaction termR IS VC is simplified to −λ 27CVC andC C1 is constant at 0.1. Figure 9 summarises the results of these two numerical experiments with the initial conditionC C1 = 0.1. From Fig. 9a, we can infer that convection marginally dominates diffusion in the centre of the embryo, but diffusion still dominates in the periphery. Figure 9b and c is very similar to their counterparts in Fig. 5 in terms of the baseline of C VC , its spatial variations, and the steepness of its gradients. When VEGFC does not bind to collagen I, its spatiotemporal dynamics in Fig. 9d are completely different. The interstitial flow pushes most of the embryo's VEGFC to the periphery. The baseline of C VC is two orders of magnitude higher than when VEGFC is sequestered by collagen I. However, the LEC is in a VEGFC gradient which only decreases twofold from the PCV to the DA. This shallow gradient is unsuitable for morphogenetic functions. Its directionality does not allow VEGFC to guide the migrating LEC to the horizontal myoseptum either.
We hypothesise that the pressure field is the key regulator of VEGFC's spatiotemporal dynamics in the last numerical experiment. To test this hypothesis, we will repeat it with a steeper and asymmetric pressure field which is given bỹ  Figure 10 illustrates the resulting dynamics ofC VC . Steepening the pressure gradient from the DA to the PCV steepens the gradient ofC VC in that region too. It is now a threefold gradient, so VEGFC may act as a morphogen for the LEC. The pressure drop from the DA to the PCV is larger than that from the DA to the DLAV. The interstitial flow will thus be faster on the ventral side of the DA. VEGFC goes with the interstitial flow, forming an asymmetric concentration field. If VEGFC is a chemorepellent, it can guide the LEC down the gradient towards the horizontal myoseptum. Considering VEGFC promotes survival, proliferation, and migration in LECs, it is unlikely to be one. Nevertheless, the pressure field in the numerical experiment is arbitrary. Reversing its direction will allow VEGFC to chemoattract the LEC to its destination.
To conclude, there is a tension between convection and VEGFC sequestration by the ECM in an embryo. Even when convection is marginally dominant in the embryo, its interstitial flow cannot influence VEGFC when the latter is sequestered by the ECM. When VEGFC is influenced by the flow, the pressure field in the embryo is the key determinant of its spatiotemporal dynamics. Its steepness determines whether VEGFC can be a morphogen for the migrating LECs in the embryo; its directionality determines whether VEGFC can guide their migration by chemotaxis.

Channelisation
Section 3.3 concludes that the spatial variations of MMP2 are small in a diffusiondominant embryo. MMP2 will gradually render such an embryo convection-dominant without introducing significant spatial effects. The embryo described in Sect. 3.4 is convection-dominant through its initialC C1 . Section 3.4 concludes that asymmetric and steep concentration gradients are possible after the manual switch from diffusion to convection. The prerequisites are an asymmetric and steep pressure field and an ECM which does not sequester the diffusing species under consideration. However, we have not considered how MMP2 fits into the picture after it takes a diffusion-dominant embryo into its convection-dominant regime. It makes us wonder whether MMP2 will behave differently in a convection-dominant embryo. To investigate this possibility, we will perform a numerical experiment. Apart from the interstitial flow Eqs. (27) and (48), we will solve the governing equations for MMP2, VEGFC, and collagen I in this numerical experiment. MMP2 will be produced at a constant rate like the numerical experiment in Sect. 3.3. The reaction terms areR LEC M2 = 10,R IS M2 = −λ 15CM2 , . We will also use the asymmetric pressure field given by equations (56), (57), and (58). We will not model MMP2-TIMP2 binding because MMP2 and TIMP2 are produced at the same source and have similar diffusion coefficients, so this phenomenon will only affect the baseline concentration of MMP2, not its spatial distribution. Figure 11 summarises the simulation results of this numerical experiment. There is a positive feedback loop between convection and ECM remodelling. Due to the asymmetric interstitial flow, MMP2 is concentrated at the ventral end of the embryo. This asymmetry inC M2 results in preferential degradation of collagen I at the ventral end, resulting in a gradient ofC C1 which decreases ventrally from the DA. This gradient will strengthen the interstitial flow in that direction to complete the feedback loop. The impact of this loop on VEGFC is shown in Fig. 11d. Its concentration gradient steepens with time because it also follows the ever increasing interstitial flow. If VEGFC is a morphogen and chemotactic factor for LECs, this positive feedback loop enhances these functions.
Our idealised zebrafish embryo is too thin for any spatial effects to show up in the x-direction. Considering our dimensions are only estimates, we will repeat the last numerical experiment with a triply widened embryo. The modified geometry will go fromx = −0.15 (outer boundary) tox = 0 (line of symmetry). The simulation results are in Fig. 12. Figure 12a, between the DA and the PCV, clearly shows some variations inC C1 in the x-direction. Therefore, we will focus on the cut line fromx = −0.15 tõ x = 0 atỹ = −0.23, a point between the LEC and the PCV. Figure 12b illustrates that a channel with lowC C1 forms gradually in the plane where the LEC and the blood vessels are located. Figure 12c indicates that VEGFC accumulates in this channel. We know that VEGFC is a growth factor for LECs, while collagen I is a physical barrier to their migration. Therefore, they are likely to stay in the channel. Channelisation ensures the LECs will migrate in the y-direction and the lymphatic vessels they form will lie in the same plane as the blood vessels.
To conclude, an asymmetric interstitial flow and the collagenolytic action of MMP2 form a positive feedback loop in a zebrafish embryo. The outcome is a channel of low C C1 and highC VC near the blood vessels' plane. The channel can keep the embryo's LECs migrating on a ventral-dorsal path, so the resulting lymphatic vessels will stay close to the blood vessels.

Concentration Gradients Vanish When Collagen I is Insufficient
Finally, we will consider the case where the initialC C1 is 0.01. Since there is so little collagen I to begin with, we will not consider its degradation by MMP2 and with- out MMP2, and we will not consider TIMP2 either. Therefore, the modified model will only include the interstitial flow Eqs. (27) and (48) and the governing equations for VEGFC, collagen I, and VEGFC bound to collagen I. To control for the kinetic properties in our model, we will increase the binding terms of VEGFC and collagen I by two orders of magnitude. Therefore, the reaction terms areR IS We can further simplify the effective diffusivity of VEGFC too. The term inside the brackets in λ 1,VC 1 − λ 2,VC λ 3CC1 is between 0.94 and 1 now. Approximating it as unity, we will retain λ 1,VC only. Figure 13a and b summarises the simulation results of this numerical experiment. The |ũ|Pe profile shows that convection is orders of magnitude more important than diffusion in the central region of the embryo. As a result, the interstitial flow flushes all the VEGFC in the embryo to the periphery where it forms boundary layers. These The spatiotemporal dynamics of VEGFC on the cut line running fromỹ = −0.45 toỹ = 0.45 atx = 0. In this numerical experiment, VEGFC is sequestered by collagen I. c The same dynamics as (b) when VEGFC does not interact with collagen I gradients cannot influence the LEC in the central region. We will repeat the numerical experiment by neglecting the interactions between VEGFC and collagen I. This means only the interstitial flow equations and the governing equation for VEGFC will be solved,R IS VC = −λ 27CVC , andC C1 is constant at 0.01. Figure 13c shows the resulting spatiotemporal dynamics of VEGFC. The central regions of Fig. 13b and c are the same, proving that sequestration by collagen I cannot protect VEGFC from such a powerful interstitial flow.
To conclude, when collagen I is insufficient, the interstitial flow in a zebrafish embryo will clear its central region of VEGFC. As a result, VEGFC cannot form a suitable concentration gradient to guide the embryo's LECs as a morphogen or a chemotactic factor. When the interstitial flow is strong enough, sequestration by collagen I cannot protect VEGFC from going with the flow. However, our pressure field is only an estimate. The ECM contains components other than collagen I too.
Our hydraulic conductivity and interstitial flow rate in this subsection are probably overestimates and not representative of zebrafish embryos.

Discussion
Our study proposes a model of the spatiotemporal dynamics of VEGFC in the trunk of a zebrafish embryo. The model parameters come from different studies. Many of them do not have in vivo sources, and many more do not pertain to zebrafish embryos. Some of them are simply unavailable, such as our use of a reported VEGF production rate as a surrogate for the required VEGFC production rate. To complicate matters, the conditions in an embryo are dynamic, so the parameters should not be static as they are in our model. Therefore, we cannot use the model to simulate what actually happens in a zebrafish embryo. Nonetheless, our simulations are scenarios within the realms of possibility. There are three scenarios in which VEGFC can be a morphogen or chemotactic factor for the LECs migrating from the PCV to the horizontal myoseptum of a zebrafish embryo.

Scenario 1
This scenario is a diffusion-dominant embryo where VEGFC is sequestered by the ECM. Under these circumstances, VEGFC will have a concentration profile which peaks at and decreases symmetrically from its source, the DA. The gradient between the PCV and the DA is perfect for VEGFC to act as a morphogen for the migrating LECs. However, the symmetry of the profile makes VEGFC an unlikely candidate for their chemotactic factor.
In general, the combination of diffusion and sequestration by an ECM creates shortrange, steep, and symmetric gradients of mobile species.

Scenario 2
Scenario 2 is an embryo marginally dominated by convection. When VEGFC binds to the embryo's ECM, it is protected from the interstitial flow therein. VEGFC will behave like it does in scenario 1.

Scenario 3
Scenario 3 is also an embryo marginally dominated by convection, but VEGFC does not bind to the ECM. In this case, the pressure field is the key factor regulating the behaviour of VEGFC. A sufficiently steep pressure field will establish a VEGFC concentration gradient appropriate for morphogenetic functions. An asymmetric pressure field will give a directionality to the gradient, thus allowing VEGFC to guide cell migration by chemotaxis.
In general, a steep and asymmetric pressure field will create steep, asymmetric, and embryo-wide concentration gradients of mobile species. The prerequisites are convection being marginally dominant and a nonbinding ECM.

Collagen I Degradation is a Control Mechanism
Collagen I concentration is the key variable that determines the dominant mode of transport, hence the scenario which applies to an embryo. Collagen I degradation by MMP2 is therefore an important means to control the behaviour of VEGFC. There are two ways MMP2 can exert its influence.
First, it allows an embryo to switch back and forth between the three scenarios. When diffusion is dominant, the distribution of MMP2 in the embryo is almost homogeneous. Collagen I will degrade almost uniformly until convection takes over. Conversely, although our model cannot demonstrate this feature, the embryo can lower MMP2 production and increase collagen I production, thus tilting the balance in diffusion's favour. This control mechanism is useful for patterning the embryo. The embryo can start in scenario 3 where embryo-wide morphogen gradients divide it into segments of different cell types. Then, it can switch to scenario 1 where local gradients will fine-tune the development of each segment.
Second, MMP2 and an asymmetric interstitial flow can form a positive feedback mechanism in scenario 3. The flow goes asymmetrically on the ventral-dorsal axis. It will concentrate MMP2 in one direction on this axis. Collagen I will as a result degrade preferentially in this direction, further skewing the interstitial flow. VEGFC will also go with the flow, thus steepening its own concentration gradient. The increased steepness makes VEGFC a more effective morphogen, while the increased asymmetry makes it a more effective chemotactic factor. Because collagen I degrades preferentially on one axis, a channel of abundant VEGFC and scarce collagen I will form along this axis. Because VEGFC is a growth factor and collagen I is a barrier to cell migration, the embryo's LECs will migrate in the channel. Channelisation therefore ensures LEC migration occurs on the ventral-dorsal axis.
Finally, when collagen I is insufficient, the embryo will be overwhelmingly dominated by convection. The interstitial flow will be so strong that no gradients of mobile species can form in the embryo. This scenario is unlikely to occur in an embryo, however.

Future Work
The developmental processes of living organisms are shaped by evolution. The fact that VEGFC is distributed like a morphogen in multiple scenarios suggests that it is actually one. On the other hand, even if VEGFC is a chemoattractant or chemorepellent, it can only perform its function in an asymmetric pressure field in scenario 3. A different molecular species is probably required to guide the PCV-derived LECs to the horizontal myoseptum in a zebrafish embryo. However, to validate our predictions, we need to comprehend the intracellular responses of an LEC to a concentration gradient of VEGFC.
Even if a species other than VEGFC induces the PCV-derived LECs to differentiate, it probably forms an appropriate concentration profile in one of the scenarios this study proposes. The proposed control mechanisms probably regulate its spatiotemporal dynamics too. More generally, these scenarios and mechanisms are frameworks within which embryonic development and tissue regeneration can be understood.
The mathematical model per se is a valuable tool for mathematical and theoretical biologists. To the best of our knowledge, it is the first mathematical model that considers an interstitial flow, a remodelling ECM, and the ECM's ability to sequester a mobile species together. Its solution costs very little computationally, so it can be extended to more complex geometries, larger domains, and higher dimensions. By changing the model geometry, its biochemical reaction network, and the model parameters, one can simulate the dynamics of a diverse range of biological phenomena.