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.

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

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 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 (Cha et al. 2012). Further guidance cues for the LECs are thought to be provided by the motor neuron axons positioned along the horizontal myoseptum (Cha et al. 2012). 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 (Cha et al. 2012).

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 (Cha et al. 2012). 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 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.

Fig. 2

Idealised geometry of a ventral–dorsal slice of a zebrafish trunk between 36 and 48 h post-fertilisation. This figure shows the idealised geometry in the anterior–posterior view. This slice is one of the 30 slices with secondary sprouts from the posterior cardinal vein. The empty circles are, from top to bottom, the dorsal longitudinal anastomotic vessel (DLAV), the dorsal aorta (DA), and the posterior cardinal vein (PCV). The solid circle is a lymphatic endothelial cell (LEC) which has exited the posterior cardinal vein; it is halfway between the dorsal aorta and the posterior cardinal vein. The dot in the middle of the figure indicates the horizontal myoseptum, which is the destination of the LEC. In this study, we consider the LEC to be stationary

Table 1

Dimensions of the idealised geometry and its internal structures

Quantity measured	Time	Measurement (μm)	References
Total height	96 HPF	434	McGee et al. (2012)
Total width	72 HPF	43	Hermans et al. (2010)
PCV diameter	96 HPF	20	Coffindaffer-Wilson et al. (2011b)
DA diameter	96 HPF	12	Coffindaffer-Wilson et al. (2011b)
DLAV diameter	96 HPF	13	Coffindaffer-Wilson et al. (2011b)
PCV-DA distance	96 HPF	51	Coffindaffer-Wilson et al. (2011b)
DA-DLAV distance	96 HPF	151	Coffindaffer-Wilson et al. (2011b)
LEC diameter	–	10	Yaniv et al. (2006)

PCV posterior cardinal vein, DA dorsal aorta, DLAV dorsal longitudinal anastomotic vessel, LEC lymphatic endothelial cell, HPF hours post-fertilisation

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 zebrafish begins by 30 HPF (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, () the specific hydraulic conductivity of the ECM, and (µm/s) the interstitial fluid velocity, we can write Brinkman’s equation asThere are two dependent variables in Eq. (1): P and , 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 byTo 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 as the position vector in our geometry and 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 () to represent the density of the interstitial fluid, (cm/Pa/s) the DA vascular permeability, and (mmHg) the pressure inside the DA. Mathematically, the relation gives the mass flux from the DA surface () into the interstitial space asWe can derive the boundary conditions on the PCV and DLAV surfaces ( and ) 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 , and the LEC surface, seen from the interstitial space domain, . 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, 1997; Siegfried 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 MMP2-related 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.

Fig. 3

Biochemical reaction network of the model. M2P, proMMP2; M2, MMP2; T2; TIMP2; C1, collagen I; MT1, MT1-MMP. A dot between two species means they are complexed together in one molecule. Only proMMP2, MMP2, and TIMP2 are present in both domains and can cross the boundary between them. A mobile species undergoes diffusion and/or convection; an immobile one does not. Only the red events are represented by the mathematical model developed in this paper

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 (M) to represent the molar concentration of species i; t (s), time; (µ/s), the effective diffusivity of species i; , the volume fraction where diffusion occurs; (µm/s), the velocity from our mathematical model of the interstitial flow; (M/s), the net rate of production of species i at a point in the interstitial space. The equation isDiffusion 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 () and the combined mass concentration of free and VEGFC-bound collagen I by (kg/), 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 (µ/s), the volume fraction of dry collagen I fibrils (without associated water molecules) by , the radius of a collagen I fibril by (µm), and the Stokes-Einstein radius of species i by (µm), the equation isUsing () 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 () and temperature by T (K), the Stokes–Einstein radius is given in Einstein (1905) asWe 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 (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 formIn the LEC domain again, the governing equations for the immobile species are a set of ordinary differential equations as follows,The reaction terms, and , 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, MMP2collagen 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 and Popel (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, , are given in Table 2. The reaction terms in the LEC, , are given in Table 3.

Table 2

Reaction terms in the interstitial space

Reaction term	Form	Equations
RM2PIS	-kM2PdegC_M2P	(7)
RM2IS	-kM2,T2onC_M2C_T2+kM2,T2offC_M2·T2-kM2degC_M2	(7)
RT2IS	-kM2,T2onC_M2C_T2+kM2,T2offC_M2·T2-kT2degC_T2	(7)
RM2·T2IS	kM2,T2onC_M2C_T2-kM2,T2offC_M2·T2-kM2·T2degC_M2·T2	(7)
RVCIS	-kVC,C1onC_VCC_C1+kVC,C1offC_VC·C1-kVCdegC_VC	(7)
RC1IS	-kM2,C1catC_M2C_C1KMM2,C1+C_C1-kVC,C1onC_VCC_C1+kVC,C1offC_VC·C1	(12)
RVC·C1IS	kVC,C1onC_VCC_C1-kVC,C1offC_VC·C1	(12)

M2 and T2 complex reversibly. VEGFC binds to C1 reversibly. M2 degrades C1 catalytically. M2P, M2, T2, , and VEGFC degrade in the interstitial space. () and () are the binding and unbinding rate constants of VC and C1; () and (), the binding and unbinding rate constants of M2 and T2; (), the turnover number in the degradation of C1 by M2; (M), the Michaelis–Menten constant in the degradation of C1 by M2; (), the degradation rate constant of species i. M2P, proMMP2; M2, MMP2; T2, TIMP2; VC, VEGFC; C1, collagen I

Table 3

Reaction terms in the lymphatic endothelial cell

Reaction term	Form	Equations
RM2PLEC	P_M2P-kMT1·T2,M2PonC_MT1·T2C_M2P	(13)
	+ kMT1·T2,M2PoffC_MT1·T2·M2P
RM2LEC	kacteffC_MT1·T2·M2PC_MT1	(13)
RT2LEC	P_T2-kMT1,T2onC_MT1C_T2+kMT1,T2offC_MT1·T2	(13)
RMT1LEC	-kMT1,T2onC_MT1C_T2+kMT1,T2offC_MT1·T2	(14)
RMT1·T2LEC	kMT1,T2onC_MT1C_T2-kMT1,T2offC_MT1·T2	(14)
	-kMT1·T2,M2PonC_MT1·T2C_M2P
	+ kMT1·T2,M2PoffC_MT1·T2·M2P
	+ kacteffC_MT1·T2·M2PC_MT1
RMT1·T2·M2PLEC	kMT1·T2,M2PonC_MT1·T2C_M2P	(14)
	-kMT1·T2,M2PoffC_MT1·T2·M2P
	-kacteffC_MT1·T2·M2PC_MT1

M2P and T2 are produced at constant rates in the lymphatic endothelial cell. T2 binds to MT1 reversibly. M2P binds to MT1T2 reversibly. MT1 activates the M2P in MT1T2M2P to form M2 and release MT1 T2. and () are the production rates of M2P and T2 by the lymphatic endothelial cell; () and (), the binding and unbinding rate constants of MT1 and T2; () and (), the binding and unbinding rate constants of and M2P; (), the activation rate constant of M2. M2P, proMMP2; M2, MMP2; T2, TIMP2; MT1, MT1-MMP

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 four outer boundaries, the no-flux boundary conditions for proMMP2, MMP2, TIMP2, MMP2TIMP2, and VEGFC are given byOn the PCV and the DLAV, they are On the DA, the boundary conditions for the mobile species except VEGFC areFor VEGFC, with (mol µ) being the release rate of VEGFC on the surface of the DA, the constant flux isOn 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 . Seen from the LEC domain, the surface is . Continuity applies to proMMP2, MMP2, and TIMP2 as follows, On the LEC surface seen from the interstitial space, we will apply no-flux conditions on MMP2TIMP2 and VEGFC,We also need a set of initial concentrations. In the LEC domain, we will assume that only MT1-MMP is present at ; in the interstitial space, we will assume that only collagen I is present at . We will label the initial concentrations of these species by and , 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 () as the hydraulic conductivity of the ECM and 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 byWith () being the molar mass of collagen I, is given byEquation (24) assumes a density of 1 kg 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 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 .

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.

Table 4

Parameters of the interstitial flow component of the mathematical model

Parameter	Definition	Value	Reference
P^DA	DA pressure	0.1844 mmHg	Hu et al. (2000)
P^DLAV	DLAV pressure	0 mmHg	Assumed
P^PCV	PCV pressure	0 mmHg	Assumed
L_DA	DA permeability	7.20×10^6cms^-1cmH_2O^-1 cm s^-1 cmH_2O^-1	Jain (1987)
L_DLAV	DLAV permeability	7.20×10^6cms^-1cmH_2O^-1	Jain (1987)
L_PCV	PCV permeability	7.20×10^6 cm s^-1 cmH_2O^-1	Jain (1987)
M_C1	Molar mass of collagen I	300 kg mol^-1	Karagiannis and Popel (2006)
μ	IF dynamic viscosity	1.200 cP	Swartz and Fleury (2007)
ρ	IF density	1025 kg m^-3	Frcitas (1998)

DA dorsal aorta, PCV posterior cardinal vein, DLAV dorsal longitudinal anastomotic vessel, IF interstitial fluid

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.

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, cm cm. 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 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 .

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 , the dynamic viscosity of blood plasma is 1.2 cP (Swartz and Fleury 2007) and its density is 1025 (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 is 2 nm (Karagiannis and Popel 2006). We can also find the required 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, and , are 0.75 and 1.89 (Levick 1987). These parameters are summarised in Table 5.

Table 5

Transport parameters in the reaction–diffusion–convection equation and its simplified forms

Parameter	Definition	Value	References
DVC∞	Diffusivity of VC	5.01×10^-7 cm^2/s	Berk et al. (1993)
DM2∞	Diffusivity of M2	0.85×10^-6 cm^2/s	Karagiannis and Popel (2006)
DM2P∞	Diffusivity of M2P	0.80×10^-6 cm^2/s	Karagiannis and Popel (2006)
DM2·T2∞	Diffusivity of M2·T2	0.75×10^-6 cm^2/s	Karagiannis and Popel (2006)
DT2∞	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)
v_C1h	Specific volume of hydrated C1	1.89 cm^3/s	Levick (1987)
T	Temperature	298 K	Assumed

M2P, proMMP2; M2, MMP2; T2, TIMP2; VC, VEGFC; C1, collagen I

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 using the known LEC diameter of 10 . 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.

Table 6

Kinetic parameters in the reaction–diffusion–convection equation and its simplified forms

Parameter	Value	References
kVC,C1on	3.60×10^4 M^-1s^-1	Köhn-Luque et al. (2013)
kVC,C1off	3.60×10^-3 s^-1	Köhn-Luque et al. (2013)
kM2,T2on	5.90×10^6 M^-1 s^-1	Karagiannis and Popel (2004)
kM2,T2off	6.30 s^-1	Karagiannis and Popel (2004)
kMT1,T2on	3.54×10^6 M^-1 s^-1	Toth et al. (2002)
kMT1,T2off	2×10^-4 s^-1	Toth et al. (2002)
kMT1·T2,M2Pon	0.14×10^6 M^-1 s^-1	Karagiannis and Popel (2004)
kMT1·T2,M2Poff	4.70×10^-3 s^-1	Karagiannis and Popel (2004)
kacteff	2.80×10^3 M^-1 s^-1	Karagiannis and Popel (2004)
kM2,C1cat	4.50×10^-3 s^-1	Karagiannis and Popel (2004)
KMM2,C1	8.50×10^-6 M	Karagiannis and Popel (2004)
kVCdeg	10^-4 s^-1	Hashambhoy et al. (2011)
kM2deg	10^-4 s^-1	Hashambhoy et al. (2011)
kM2Pdeg	10^-4 s^-1	Hashambhoy et al. (2011)
kM2·T2deg	10^-4 s^-1	Hashambhoy et al. (2011)
kT2deg	10^-4 s^-1	Hashambhoy et al. (2011)
P_M2P	2.64×10^-8 M s^-1	Vempati et al. (2010)
P_T2	1.54×10^-10 M s^-1	Vempati et al. (2010)
RVCDA	1.65×10^-17 mol dm^-2 s^-1	Hashambhoy et al. (2011)

means the binding rate constant of species i and j; , their unbinding rate constant; , the activation rate constant of M2; , the turnover number in the degradation of C1 by M2; , the Michaelis–Menten constant in the degradation of C1 by M2; , the degradation rate constant of species i; , the production rate of species i; , the production rate of VC on the surface of the dorsal aorta

M2P, proMMP2; M2, MMP2; T2, TIMP2; VC, VEGFC; C1, collagen I; MT1, MT1-MMP

Finally, we need the initial concentration of MT1-MMP in the LEC and that of collagen I in the interstitial space. is unavailable for endothelial cells. However, a value based on other cell types, 180,000 molecules/cell or M, is in Karagiannis and Popel (2006). In adult tissues, the concentration of collagen I ranges from to M (Levick 1987; Karagiannis and Popel 2006). We will use the midpoint of this range, so is M for adult tissues. The collagen content of frog embryos and larvae ranges from M to M (Edds Jr 1958). We will use the midpoint of this range, 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, , we need the length scale, L; we will use the largest dimension of the geometry, 434 µm. We will use mmHg to scale the pressure field, . Because we are interested in the period from 36 to 48 HPF, we will use a time scale, , of 12 h. The nondimensionalised time is therefore . Velocity is scaled like . The concentrations are nondimensionalised like . 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, . We will use the adult value, M. The concentrations of MT1-MMP and its two complexes will always add up to the initial concentration of MT1-MMP, M. We will use this value for , , and . 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.

Table 7

Scales used for nondimensionalisation

Scale	Description	Value
C_C1,s	Concentration scale for C1	3.50×10^-4 M
C_VC,s	Concentration scale for VC	1.64×10^-10 M
C_VC·C1,s	Concentration scale for VC·C1	8.93×10^-5 M
C_M2,s	Concentration scale for M2	3.94×10^-5 M
C_M2P,s	Concentration scale for M2P	1.14×10^-3 M
C_M2·T2,s	Concentration scale for M2·T2	6.68×10^1 M
C_MT1,s	Concentration scale for MT1	5.71×10^-7 M
C_MT1·T2,s	Concentration scale for MT1·T2	5.71×10^-7 M
C_MT1·T2·M2P,s	Concentration scale for MT1·T2·M2P	5.71×10^-7 M
C_T2,s	Concentration scale for T2	6.65×10^-6 M
L	Length scale	434 μm
P^DA	Pressure scale and DA pressure	0.1844 mmHg
τ	Time scale	43200 s
U	Velocity scale	1.371×10^-4 μm/s

M2P, proMMP2; M2, MMP2; T2, TIMP2; VC, VEGFC; C1, collagen 1; MT1, MT1-MMP

Table 8

Dimensionless parameters in the nondimensionalised interstitial flow and reaction–diffusion–convection equation

Parameter	Form	Value
α	Constant	−2.70
η_1	μULP^DAβ(M_C1C_C1,s)^α	1
η_2	C_VC·C1,sC_C1,s	0.255
η_3	μULP^DA	1.542×10^-11
η_DA	L_DAP^DAU	1.317×10^14
η_PCV	L_PCVP^DAU	1.317×10^14
η_DLAV	2L_DLAVP^DAU	2.634×10^14
λ_1,VC	DVC∞τL^2	1.15×10^1
λ_1,M2	DM2∞τL^2	1.95×10^1
λ_1,M2P	DM2P∞τL^2	1.83×10^1
λ_1,M2·T2	DM2∞·T2τL^2	1.72×10^1
λ_1,T2	DT2∞τL^2	2.52×10^1
λ_2,VC	k_BT6πμDVC∞r_f	1.81
λ_2,M2	k_BT6πμDM2∞r_f	1.07
λ_2,M2P	k_BT6πμDM2P∞r_f	1.14
λ_2,M2·T2	k_BT6πμDM2∞·T2r_f	1.21
λ_2,T2	k_BT6πμDT2∞r_f	8.26×10^-1
λ_3	v_C1M_C1C_C1,s	7.88×10^-2
λ_4	v_C1M_C1C_VC·C1,s	2.01×10^-2
λ_5	v_C1hM_C1C_C1,s	1.98×10^-1
λ_6	v_C1hM_C1C_VC·C1,s	5.06×10^-2
λ_7	UτL	1.36×10^-2

DA dorsal aorta, PCV posterior cardinal vein, DLAV dorsal longitudinal anastomotic vessel

M2P, proMMP2; M2, MMP2; T2, TIMP2; VC, VEGFC; C1, collagen I

Table 9

Dimensionless parameters in the nondimensionalised reaction terms

Parameter	Form	Value
λ_8	P_M2PτC_M2P,s	1
λ_9	kMT1·T2,M2PonτC_MT1·T2,s	3.45×10^3
λ_10	kMT1·T2,M2PoffτC_MT1·T2·M2P,sC_M2P,s	1.02×10^-1
λ_11	kM2Pdegτ	4.32
λ_12	kacteffτC_MT1·T2·M2P,sC_MT1,sC_M2,s	1
λ_13	kM2,T2onτC_T2,s	1.69×10^6
λ_14	kM2,T2offτC_M2·T2,sC_M2,s	4.61×10^11
λ_15	kM2degτ	4.32
λ_16	P_T2τC_T2,s	1
λ_17	kMT1,T2onτC_MT1,s	8.73×10^4
λ_18	kMT1,T2offτC_MT1·T2,sC_T2,s	7.42×10^-1
λ_19	kM2,T2onτC_M2,s	1×10^7
λ_20	kM2,T2offτC_M2·T2,sC_T2,s	2.73×10^12
λ_21	kT2degτ	4.32
λ_22	kM2,T2onτC_M2,sC_T2,sC_M2·T2,s	1
λ_23	kM2,T2offτ	2.72×10^5
λ_24	kM2·T2degτ	4.32
λ_25	kVC,C1onτC_C1,s	5.44×10^5
λ_26	kVC,C1offτC_VC·C1,sC_VC,s	8.47×10^7
λ_27	kVCdegτ	4.32
λVCDA	RVCDAτC_VC,sL	1
λ_28	kM2,C1catτC_M2,s	7.66×10^-3
λ_29	kVC,C1onτC_VC,s	2.55×10^-1
λ_30	kVC,C1offτC_VC·C1,sC_C1,s	3.97×10^1
λ_31	kVC,C1onτC_VC,sC_C1,sC_VC·C1,s	1
λ_32	kVC,C1offτ	1.56×10^2
λ_33	kMT1,T2onτC_T2,s	1.02×10^6
λ_34	kMT1,T2offτC_MT1·T2,sC_MT1,s	8.64
λ_35	kMT1,T2onτC_MT1,sC_T2,sC_MT1·T2,s	1.02×10^6
λ_36	kMT1,T2offτ	8.64
λ_37	kMT1·T2,M2PonτC_M2P,s	6.89×10^6
λ_38	kMT1·T2,M2PoffτC_MT1·T2·M2P,sC_MT1·T2,s	2.03×10^2
λ_39	kacteffτC_MT1·T2·M2P,sC_MT1,sC_MT1·T2,s	6.91×10^1
λ_40	kMT1·T2,M2PonτC_MT1·T2,sC_M2P,sC_MT1·T2·M2P,s	6.89×10^6
λ_41	kMT1·T2,M2Poffτ	2.03×10^2
λ_42	kacteffτC_MT1,s	6.91×10^1

It should be noted that is not really dimensionless and is in M, but the term is dimensionless because the denominator is also in M

M2P, proMMP2; M2, MMP2; T2, TIMP2; VC, VEGFC; C1, collagen 1; MT1, MT1-MMP; DA dorsal aorta

First, we will nondimensionalise the interstitial flow equations. In this study, a tilde represents a nondimensionalised variable, for example, . Combining Eqs. (23) and (24), we can rearrange the resulting equation to write in terms of and . We can convert to using . Defining the constants and , we can write as . Substituting this equation for into Eq. (1) and nondimensionalising the variables, we can write the nondimensionalised Brinkman’s equation as follows,To simplify Eq. (25), we will lump the parameters into dimensionless groups, , , and . To determine the scale of , we can use either or . Choosing will lead to µm/s; choosing will lead to m/s. The velocity for an interstitial flow is reported to range from 0.1 to 2 µm/s (Swartz and Fleury 2007), so gives a more appropriate velocity scale. This choice also leads to . Therefore, Eq. (25) becomesNondimensionalising Eq. (2) and the boundary conditions, then defining the dimensionless groups , , and , we will obtain the remaining nondimensionalised interstitial flow equations and their boundary conditions, Second, we need to nondimensionalise the equations governing the concentration fields in the interstitial space, (7) and (12). To do so, we need to introduce , , and . Then, the nondimensionalised reaction–diffusion–convection equation takes the formWe will define the dimensionless parameters and . Then, we can write as can be expressed in terms of the dimensionless parameters , , , and . This expression isFor the immobile species in the interstitial space, the reaction–diffusion–convection equation is reduced toThird, we will nondimensionalise the equations governing the concentration fields in the LEC domain, (13) and (14). After introducing , we can write downFor the immobile species in the LEC domain, the reaction–diffusion equation is reduced toFourth, we will nondimensionalise the boundary and initial conditions of the concentrations. Equations (15–17), which apply to proMMP2, MMP2, TIMP2, MMP2TIMP2, and VEGFC, have the following nondimensionalised forms, The boundary conditions for proMMP2, MMP2, TIMP2, and MMP2TIMP2 on the DA’s surface, represented by Eq. (18), have the following nondimensionalised form,With being the dimensionless production rate of VEGFC, we can scale by setting , so M. Then, the influx of VEGFC from the DA surface into the interstitial space, Eq. (19), has the nondimensionalised formThe 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 MMP2TIMP2 and VEGFC.

The initial concentrations of MT1-MMP in the LEC domain and collagen I in the interstitial space domain are and , respectively.

Finally, we need to define the nondimensionalised reaction terms, and , in terms of nondimensionalised concentrations. We need the remaining concentration scales. We will scale by setting the dimensionless production term of proMMP2 to unity, . This leads to M. We will set the dimensionless production term of MMP2 to unity, . This gives the concentration scale for MMP2, M. Setting the dimensionless production term of TIMP2 to unity, , we can obtain the concentration scale M. If the binding term of MMP2 and TIMP2 is set to unity, , the concentration scale is M. Similarly, by equating the binding term of VEGFC and collagen I to unity, , we can scale the concentration of bound VEGFC, M. With all the scales determined, the reaction terms can easily be worked out. They are summarised in Table 10.

Table 10

Nondimensionalised reaction terms in the interstitial space and lymphatic endothelial cell domains

Reaction term	Form	Equations
R~M2PIS	-λ_11C~_M2P	(32)
R~M2IS	-λ_13C~_M2C~_T2+λ_14C~_M2·T2-λ_15C~_M2	(32)
R~T2IS	-λ_19C~_M2C~_T2+λ_20C~_M2·T2-λ_21C~_T2	(32)
R~M2·T2IS	C~_M2C~_T2-λ_23C~_M2·T2-λ_24C~_M2·T2	(32)
R~VCIS	-λ_25C~_VCC~_C1+λ_26C~_VC·C1-λ_27C~_VC	(32)
R~C1IS	-λ_28C~_M2C~_C1KMM2,C1+C_C1,sC~_C1-λ_29C~_VCC~_C1+λ_30C~_VC·C1	(35)
R~VC·C1IS	C~_VCC~_C1-λ_32C~_VC·C1	(35)
R~M2PLEC	1-λ_9C~_MT1·T2C~_M2P+λ_10C~_MT1·T2·M2P	(36)
R~M2LEC	C~_MT1·T2·M2PC~_MT1	(36)
R~T2LEC	1-λ_17C~_MT1C~_T2+λ_18C~_MT1·T2	(36)
R~MT1LEC	-λ_33C~_MT1C~_T2+λ_34C~_MT1·T2	(37)
R~MT1·T2LEC	λ_35C~_MT1C~_T2-λ_36C~_MT1·T2	(37)
	-λ_37C~_MT1·T2C~_M2P+λ_38C~_MT1·T2·M2P
	+λ_39C~_MT1·T2·M2PC~_MT1
R~MT1·T2·M2PLEC	λ_40C~_MT1·T2C~_M2P-λ_41C~_MT1·T2·M2P	(37)
	-λ_42C~_MT1·T2·M2PC~_MT1

M2P, proMMP2; M2, MMP2; T2, TIMP2; VC, VEGFC; C1, collagen 1; MT1, MT1-MMP

IS interstitial space, LEC lymphatic endothelial cell

Simplification

After nondimensionalisation, we can spot several opportunities for simplifying the model.

First, the nonlinear term in Eq. (26), , can be linearised if and differ by orders of magnitude. After factorisation, the nonlinear term becomes , which can be expanded as a Taylor series at . The result isThis expansion assumes that is true in the interstitial space domain. The binding rate constant of VEGFC and collagen I in is ; the unbinding rate constant of the same process is . The initial concentration of collagen I is unity, but it is zero for VEGFC bound by collagen I. It follows that must be orders of magnitude smaller than all the time. We also know that . 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 and define the simulated pressure, , resulting inSecond, the blood vessels are very leaky. In the equations representing the transvascular fluxes of interstitial fluid, (28–30), the dimensionless vascular permeabilities, , , and , 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 to Third, we can linearise the diffusive term in Eq. (32), . We will start with . We can apply around (Abramowitz and Stegun 1964) to obtainThis approximation assumes that holds in the interstitial space. Since and are true, must hold too. With respect to the second term, we know that and , so it is orders of magnitude less than unity. The assumption is a valid one. Also, the terms beyond are negligible. The result is .

The effective diffusivity can be simplified along similar lines. The Taylor series expansion of a general exponential function around is We can apply it to the exponential term in , resulting inAs usual, we need the magnitude of the exponent to be much smaller than one for Taylor expansion to work. Since and are true, the first term in the square root is much smaller than one; and mean the same for the second term. The square root is on the order of or smaller. Since is on the order of unity, this expansion is justified. For the same reason, the terms beyond in the expansion are negligible, thus reducing the diffusivity to .

The square root in can be simplified even further. Rewriting it as and knowing that , we will expand it as a Taylor series around ,Because the terms beyond unity are negligible, the effective diffusivity is simplified to .

The simplified Eq. (32) isThe simplified expressions and will also be applied to the boundary conditions.

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 ; 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 from to . 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 µm/s, it takes a particle s to transverse 1 by convection. From , the highest diffusivity in our model, we know that the equivalent time is s for diffusion. Nondimensionalising the shorter travelling time, s becomes . We will therefore set the initial time step at .

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 to 0.037 in nondimensionalised length.

Fig. 4

(Color figure online) Convergence plot for the mesh refinement study. It shows how the maximum velocity magnitude at changes with the number of mesh elements. The data points are the numerical results obtained with, from left to right, the ‘fine’, ‘finer’, ‘extra fine’, and ‘extremely fine’ mesh settings in COMSOL Multiphysics version 5.2

We will next optimise the maximum size of a time step. The equations governing , (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 minimum at is 0.9991. When the maximum time step is halved to 0.01, the minimum at is 0.99909. The difference is 0.001%. The biologically relevant range of 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 . This expression is based on the characteristic velocity scale, but the velocity magnitude varies in space. A better measure is , the maximum of which is 0.14909 at for VEGFC. Of all the diffusible species in our model, VEGFC has the smallest , so can only be smaller for another species. Diffusion is the dominant mode of transport in the primary simulation.

The spatiotemporal dynamics of are shown in Fig. 5. Throughout our time frame of interest, 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.

Fig. 5

(Color figure online) Spatiotemporal dynamics of VEGFC in the primary simulation. a shows the full concentration profiles of VEGFC at different time points. b Defines the coordinate system of the geometry and a cut line in the y-direction. The cut line, which is in red, goes from to at . c The concentration profiles of VEGFC on this cut line at different time points. The gaps are, from left to right, the posterior cardinal vein, the lymphatic endothelial cell, the dorsal aorta, and the dorsal longitudinal anastomotic vessel

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. At , 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 . It is estimated in Gurdon and Bourillot (2001) that a cell can read a threefold change in a morphogen’s concentration over 30 . 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 to M (Gurdon and Bourillot 2001). in the primary simulation is on the order of 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 tenfold in a numerical experiment, the baseline of will increase tenfold, but its profile’s shape will remain unchanged. This is illustrated in Fig. 6a.

Fig. 6

(Color figure online) Spatiotemporal dynamics of VEGFC in two numerical experiments. Both sets of results are shown on the cut line running from to at . a The simulated dynamics when the production rate of VEGFC is increased tenfold. b The dynamics simulated without VEGFC–collagen I interactions

We will perform another numerical experiment by setting , , , , , and . 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 of 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-MMPTIMP2proMMP2. 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.

Fig. 7

(Color figure online) Behaviour of MMP2 in the primary simulation. a The temporal dynamics of MT1-MMP and its complexes inside the lymphatic endothelial cell. The maximum concentration of each species in the cell is plotted at each time point. b The production rate of MMP2 in the cell, . The maximum production rate in the cell is plotted at each time point. c The concentration profiles of MMP2 at selected time points. d The spatiotemporal dynamics of collagen I on the cut line running from to at . MT1, MT1-MMP; T2, TIMP2; M2P, proMMP2

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, , , , and remain. Furthermore, we will use and in the modified model. In Fig. 7b, the maximum is around 0.2, so 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%.

Fig. 8

(Color figure online) Behaviour of MMP2 in a numerical experiment which replaces the activation mechanism of MMP2 with a constant production rate. a The concentration profiles of MMP2 at selected time points. b The spatiotemporal dynamics of collagen I on the cut line running from to at

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-MMP degrades the collagen I around the LECs, thereby triggering LEC migration. 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 is from M to M. Our choice of is M, so the initial condition falls at the upper end of the range.

We will carry out a numerical experiment where the initial 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 of in a diffusion-dominant zebrafish embryo. Our manual reduction in the initial has the same effect. A reduction in 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 lowering 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 become , , and .

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 reaction–diffusion–convection equation governing . The reaction term is simplified to and is constant at 0.1.

Figure 9 summarises the results of these two numerical experiments with the initial condition . 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 , 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 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.

Fig. 9

(Color figure online) Spatiotemporal dynamics of VEGFC in a marginally convection-dominant zebrafish embryo. The initial is 0.1 and collagen I does not degrade. a That convection is marginally dominant in the central region of the embryo, but diffusion is still dominant in the periphery. b The concentration profiles of VEGFC at selected time points. VEGFC is sequestered by collagen I in this case. c Plots the data of (b) on the cut line running from to at . d The spatiotemporal dynamics of VEGFC on the same cut line when VEGFC does not interact with collagen I

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 by Figure 10 illustrates the resulting dynamics of . Steepening the pressure gradient from the DA to the PCV steepens the gradient of 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.

Fig. 10

(Color figure online) Spatiotemporal dynamics of VEGFC in an asymmetric pressure field. The initial is 0.1, and collagen I does not degrade. Convection is marginally dominant in the central region of the embryo, while diffusion dominates in the periphery. VEGFC does not interact with collagen I. a A collection of VEGFC concentration profiles at selected time points. b The spatiotemporal dynamics of VEGFC on the cut line running from to at

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 diffusion-dominant 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 initial . 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 are , , , and . 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 in results in preferential degradation of collagen I at the ventral end, resulting in a gradient of 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.

Fig. 11

(Color figure online) Positive feedback loop between an asymmetric interstitial flow and the collagenolytic action of MMP2. Initially, is 0.1 and convection is marginally dominant in the central region of the embryo, but diffusion dominates in the periphery. MMP2 is produced at a constant rate to degrade collagen I. VEGFC is not sequestered by collagen I. The pressure field is asymmetric. a The spatiotemporal dynamics of MMP2. b The spatiotemporal dynamics of collagen I. c That convection becomes increasingly dominant in the region ventral to the DA. d The spatiotemporal dynamics of VEGFC on the cut line running from to at

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 from (outer boundary) to (line of symmetry). The simulation results are in Fig. 12. Figure 12a, between the DA and the PCV, clearly shows some variations in in the x-direction. Therefore, we will focus on the cut line from to at , a point between the LEC and the PCV. Figure 12b illustrates that a channel with low 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.

Fig. 12

(Color figure online) Channels form in collagen I due to an asymmetric interstitial flow and the collagenolytic action of MMP2. Initially, is 0.1 and convection is marginally dominant in the central region of the embryo, but diffusion dominates in the periphery. MMP2 is produced at a constant rate to degrade collagen I. VEGFC is not sequestered by collagen I. The pressure field is asymmetric. a The concentration profile of collagen I at . b The spatiotemporal dynamics of collagen I on the cut line running from (outer boundary) to (line of symmetry) at . This cut line is between the posterior cardinal vein and the lymphatic endothelial cell. c The spatiotemporal dynamics of VEGFC on the same cut line. d The aforementioned cut line from to at

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 and high 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 initial is 0.01. Since there is so little collagen I to begin with, we will not consider its degradation by MMP2 and without 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 are , , and . We can further simplify the effective diffusivity of VEGFC too. The term inside the brackets in is between 0.94 and 1 now. Approximating it as unity, we will retain only.

Figure 13a and b summarises the simulation results of this numerical experiment. The 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 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, , and 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.

Fig. 13

(Color figure online) Concentration gradients cannot form with insufficient collagen I. In these two numerical experiments, the initial is 0.01 and there is no collagen I degradation. a The central region of the embryo is overwhelmingly dominated by convection although diffusion dominates in the periphery. b The spatiotemporal dynamics of VEGFC on the cut line running from to at . In this numerical experiment, VEGFC is sequestered by collagen I. c The same dynamics as (b) when VEGFC does not interact with collagen I

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 short-range, 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.