Modelling the lymphatic system: challenges and opportunities

The lymphatic system is a vital part of the circulatory and immune systems, and plays an important role in homeostasis by controlling extracellular fluid volume and in combating infection. Nevertheless, there is a notable disparity in terms of research effort expended in relation to the treatment of lymphatic diseases in contrast to the cardiovascular system. While similarities to the cardiovascular system exist, there are considerable differences in their anatomy and physiology. This review outlines some of the challenges and opportunities for those engaged in modelling biological systems. The study of the lymphatic system is still in its infancy, the vast majority of the models presented in the literature to date having been developed since 2003. The number of distinct models and their variants are few in number, and only one effort has been made thus far to study the entire lymphatic network; elements of the lymphatic system such as the nodes, which act as pumps and reservoirs, have not been addressed by mathematical models. Clearly, more work will be necessary in combination with experimental verification in order to progress and update the knowledge on the function of the lymphatic system. As our knowledge and understanding of its function increase, new and more effective treatments of lymphatic diseases are bound to emerge.


INTRODUCTION
The lymphatic system constitutes a one-way transport system that operates in conjunction with the circulatory system. Its primary function is to transport excess interstitial fluid, from the interstitial space, back to the blood circulation, via the thoracic duct [1]. Along with the excess interstitial fluid, excessive proteins and waste are transported back to the circulation. The lymphatic system also acts as a conduit for immune cells and facilitates the immune response [2]. Lymph nodes across the network filter the interstitial flow and break down bacteria, viruses and waste. The lymphatic system is, therefore, not part of the circulatory system alone; rather, it is an integral part of the immune system. It plays an important role in the dissemination of cancer [3,4]. Its role in transplantation is also significant. Latest findings suggest that early lymphangiogenesis in a transplanted graft may be responsible for early rejection, but it is beneficial later on [5].
Failure of lymph drainage can be a consequence of infection, trauma, surgery, transplantation, medication or venous disease, or it may be congenital [3,6]. The build-up of interstitial fluid results in swelling and pain, and increases the risk of infection. This condition is known as lymphoedema and current treatments have limited success. Arm lymphoedema is very common in breast cancer patients and lower limb lymphoedema may lead to incapacity in severe cases.
Despite the importance of the lymphatic system in health and disease, it remains overlooked in terms of research, especially compared with the circulatory system [7,8]. Although tissue engineering has made some progress towards therapy in other areas, in the lymphatic system it is still in its infancy [4].
The lymphatic system is a highly complex system with highly variable structure and function between anatomical sites and between species. The lack of mathematical models developed for the study of its function presents a wealth of opportunities for researchers. Unfortunately, this is hindered by the lack of anatomical and physiological data.
The present review attempts to summarize the current knowledge on the lymphatic system with an emphasis on the modelling of its biomechanical behaviour. Sections 2 and 3 give a brief overview of the anatomy and physiology of the lymphatic system, respectively. Section 4 deals with the mathematical models developed thus far and discusses their key assumptions and limitations. Finally, §5 anticipates future research.

ANATOMY OF THE LYMPHATIC SYSTEM
The lymphatic system is composed of a network of vessels, termed lymphatics, lymph nodes and lymphoid organs. The interstitial fluid enters through the small lymphatic capillaries (also called initial or terminal lymphatics) that gradually combine to form larger diameter vessels, namely the pre-collectors, collectors, trunks and ducts. The interstitial fluid, which is termed lymph when inside the lymphatics, is pumped slowly by the contraction of the lymphatic vessels, which contain smooth muscle in their walls. Retrograde flow is prevented by a series of one-way valves. Figure 1 shows a simplified schematic diagram of the lymphatic system. The lymphatic network is asymmetrical: the righthand side of the head and thorax and the right arm drain into the right subclavian vein, whereas the lymphatic vessels of the rest of the body converge at the thoracic duct, which empties at the junction of the jugular and left subclavian veins.

Initial lymphatics
Initial lymphatic capillaries have a diameter of approximately 10 -60 mm with a wall thickness ranging from 50 to 100 nm and are blind-ended. Electron microscopy revealed that lymphatic capillaries are usually collapsed with irregular walls and do not have a continuous basal lamina. They comprise a monolayer of nonfenestrated endothelial cells [4]. They are supported by anchoring filaments (6 -10 nm in diameter) that keep them from completely collapsing when interstitial fluid pressure increases [10]. Their unique structure allows them to function as a one-way valve system that allows fluid to flow into them and then to close in order to stop back-flow as the internal pressure rises. Trzewik et al. [11] demonstrated this function experimentally, although they were not the first to postulate the existence of primary lymphatic valves. Leak & Burke [10] investigated the structure of the lymphatic capillaries, and speculated that they could function as valves, although at that time no proof was provided.

Pre-collecting lymphatics
The pre-collecting lymphatics connect the capillaries to the collecting vessels. They contain bicuspid one-way valves, but unlike the collecting vessels, where they are located at regular intervals, their distribution becomes more irregular towards the capillary vessels, and may comprise a single leaflet [6,12]. The pre-collectors contain one or more layers of smooth muscle cells within their walls, and are capable of performing spontaneous contractions. However, portions without muscle exist; in these parts, the endothelial layer in the pre-collectors is similar to that of initial lymphatics with a discontinuous basal lamina. Thus, these structures absorb fluid instead and therefore the precollectors have a dual role: absorption and propulsion of lymph [12].

Collecting lymphatics
The larger collecting lymphatics differ from the capillaries and pre-collectors, as they have a complete basal lamina and therefore the primary lymphatic valves are absent. They contain another type of valve that prevents retrograde flow: the secondary valve. The wall structure is similar to that of blood vessels. Three layers can be identified-the intima, media and adventitia [13]. 1 They are composed of endothelial cells, smooth muscle and collagen fibres, respectively. It has been observed that the muscle bundles are arranged in a hellicoidal manner in both collecting and pre-collecting vessels [12]. The part of the vessel between two valves is known as a lymphangion (figure 2).    [9]. 1 Lymphangions are innervated with sympathetic and parasympathetic nerves, and can perform rhythmic contractions. In humans, the lymphangions present within the head and neck have an average diameter of 0.2 mm and a length of 2 mm, yielding a diameter to length ratio of 0.1 [17]. Collecting lymphatics are approximately 1-2 mm in diameter in human lower extremities [18]. The largest (and deeper) lymphatic vessels, also referred to as trunks and ducts, have a diameter of the order of 2 mm, and the diameter to length ratio can be close to unity (figure 3). This may have implications for the modelling of the lymph flow, which will be discussed in §2.4.

Secondary valves
Generally, the lymphatic secondary valves are of the bicuspid type, regardless of species studied [19,20]. However, up to five leaflet valves as well as single leaflet ones have been observed [21]. Data regarding the size, the geometry and the number of valves are scarce, and, in Figure 2. Lymphangion. This schematic diagram cannot encompass the wide variation of vessel aspect ratio observed at different anatomical sites; and while the distension shown in the figure is perhaps exaggerated, lymphangions have been reported to double or triple their diameter during diastole [15,16]. Reproduced with permission from Quere [9].
lymphatic trunk collecting lymph vessel 2 mm (interval of the valves) f2 mm f0.2 mm Figure 3. Schematic diagram of a collecting vessel and a trunk/ duct of the human head and neck from a cadaveric study. It appears that, at this anatomical site in humans, the intervalvular distance is roughly the same. However, this is not a general case [18]. Reproduced with permission from Pan et al. [17].  humans especially, there is a lack of data on the structure and operation of the lymphatic valves [22], although valve size varies with the vessel calibre [17]. Studies of the rat spinotrapezius muscle by Mazzoni et al. [20] have shown that smooth muscle cells are absent from the valve leaflets, indicating that these are passive structures. 2 The valves comprise a monolayer of endothelial cells on a collagen matrix. Mazzoni et al. [20] suggested that the valve operation is determined by the pressure and viscous forces associated with the lymph flow. Based on the anatomy of valves studied and the presence of a buttress structure, the authors speculated that tension from the surrounding tissues may also contribute to the valve function. Recently, Davis et al. [23] investigated the valve gating in rat mesenteric lymphatics. It was observed that opening and closing of valves were occurring only because of a pressure gradient; however, the vessel muscle tone influenced the required pressure gradient. Studies of the thoracic duct of monkeys by Lee et al. [19] have shown the presence of bicuspid valves, but no buttress structure, indicating that there may be differences in the details of valve operation depending on species and anatomical site.

Lymphoid organs
Lymphoid organs are classified as being either primary or secondary organs. The primary organs are the thymus and bone marrow, and are responsible for the production and maturation of lymphocytes. The secondary organs include the spleen, Peyer's patches, appendix, the tonsils and the nodes. They are responsible for further maturation of lymphocytes and initiation of an immune response. Only the structure and function of the lymph nodes will be discussed in the context of this paper, as they are the only lymphoid organs that play a role in the active transport of lymph [9].
The structure of the node is highly complex and not easily visualized, owing to the large number of cells that are resident within it [24]. A representation of the lymph node structure is shown in figure 4a, and a more detailed cross section is shown in figure 4b.
The exterior walls of the node are covered with lymphatic smooth muscle, thus exhibiting contractile behaviour like the lymphatic vessels, although the frequency of contraction is lower. Hughes & Allen [26] and McGeown & Gallagher [27] report a frequency of 0.78 b.p.m. in bovine nodes, while Thornbury et al. [28] observed contractions at a frequency of 0.5-1 b.p.m. in ovine nodes. Tumer et al. [29] performed a frequency analysis and found the dominant rhythms to lie within the ranges of 0.01-0.04, 0.05-0.07 and 0.09-0.14 Hz. In contrast, lymphatic collecting vessels contract at a frequency in the range of 0.6-10 b.p.m. [30][31][32][33], although up to 30 b.p.m. have been reported in the literature [34]. Limited information exists on human nodal contractility, its presence and its regulating mechanisms [22].
Lymph nodes present a relatively high resistance to flow, which depends on the flow rate of lymph [35,36].
Papp et al. [36] reported that the resistance of the node in dogs was 100 times larger than that of the thoracic duct, while Browse et al. [35] found that the resistance of the nodes fell by up to six times as the flow rate in the afferent lymphatics increases. They concluded that the reduction in resistance was due to swelling of the nodes, which created more space for lymph to flow, thereby reducing their resistance. Nevertheless, despite this reduction, the nodes still impose a large resistance to the flow.
It is not surprising that the nodes have a large resistance, if one considers their internal structure. Lymph flows into the subcapsular sinus and is immediately forced to flow around the lymphoid compartment 3 and into a network of sinuses (yellow arrows in figure 4b), where it is filtered by macrophages that reside on the outer surface of the lymphoid compartment. Most of the lymph exits the node from the efferent lymphatics without entering the lymphoid compartment. While in figure 4 multiple afferents and one efferent vessel are shown, any possible combination of the number of these vessels has been observed in humans [17] as well as in mice [38]. Kowala & Schoefl [38] also observed that there is no constant relationship in the number of afferent vessels and the number of nodal compartments. Whether this variation in numbers of vessels is for overall node resistance reduction, increased lymph-filtering efficiency or both is uncertain.
The outer surface of the lymphoid compartment is covered by sinus-lining cells that create an almost impermeable membrane [25]. Large molecules cannot penetrate this layer, thus creating an effective barrier that prevents pathogens from reaching the bloodstream. Lymphocytes, lymph and smaller molecules can enter through specialized channels [24,25]. The exclusion of large molecules from the lymphoid compartment explains the well-established fact that the protein composition of efferent lymph is higher than that of the afferent one, although the exact mechanisms of mass exchange are not fully understood.
The amount of lymph entering is low and therefore the inside of the lymphoid compartment is relatively fluidfree, although controversy in the literature does exist as, according to Tretbar et al. [6], up to 50 per cent of the lymph is lost inside this compartment. Lymph flows within specialized channels, referred to as the conduit system [25,39], that start from the sinus-lining cells and end on the lumen of venules that penetrate inside the lymphoid compartment. The cells lining these venules are specially adapted to allow lymphocytes to freely enter the node from the blood and are named high endothelial venules (HEVs). HEVs also allow lymph to return to the bloodstream.

PHYSIOLOGY OF THE LYMPHATIC SYSTEM
The interstitial fluid is taken up by the small lymphatic capillaries in a process known as lymph formation. Lymph is then pumped against a pressure gradient towards the jugular vein, by rhythmic contractions of the lymphangions and by external motion of skeletal muscle, arteries and veins. This process is referred to as lymph propulsion. Lymph formation depends on external (or extrinsic), but local, effects such as interstitial fluid pressure and extracellular matrix strain and is thus a local process [3]. On the other hand, lymph propulsion is a systemic process regulated by a large number of factors: streamwise pressure gradients, transmural pressure, nerves and hormones.

Lymph formation
Several mechanisms have been proposed for the formation of lymph in initial lymphatics. Of those, only the interstitial pressure and volume theory have been widely accepted [40]. Yet, the exact regulatory mechanisms are unclear [41]. Lymph enters through openings in the endothelial layer of the capillaries, which act as one-way valves. They open when interstitial pressure increases and permit interstitial fluid flow inside the initial lymphatics. As the pressure inside the vessel rises, the valve closes. Flow inside the initial lymphatics is thought to be facilitated by fluctuations in interstitial fluid pressure and by the suction of the collecting lymphatics downstream [42]. Studies in a number of mammalian species (rodents, canines, ovines, bovines, humans) have shown no evidence of smooth muscle cells in the initial lymphatics. It follows that the driving force of lymph formation must be the result of extrinsic factors. The sole exception to this is the bat [34,40].
Lymph formation is organ dependent because of the differences in the structure and mechanical properties between tissues [43].

Lymph propulsion
After the filling of the initial lymphatics, lymph has to be pumped against a pressure gradient. Lymph propulsion is performed by the rhythmic contractions of the smooth muscles of lymphangions and facilitated by the presence of the one-way valves that prevent retrograde flow. The contractions propagate from one lymphangion to the next like a wave that causes both radial contraction and axial shortening of the vessels [33,44]. It has been reported that lymphangions rotate in vitro as well [45], most probably owing to the helical disposition of these cells within the vessel wall [12]. Although there is evidence to support that the contraction wave is transmitted between lymphangions through endothelial gap junctions [33], it has been observed that the propagation of the wave may become discontinuous. This suggests that lymphangions have a self-regulating mechanism mediated by pacemaker cells that are thought to reside downstream of the valves [31]. The contraction wave can travel in an orthograde or retrograde direction. Evidence shows that either propagation direction is equally likely to take place [33] and the volume of fluid pumped is not significantly affected by it [31].
However, there is some controversy regarding the nature of individual lymphangion contraction. The accepted view is that lymphangions contract in a peristaltic manner, that is, the contraction is a radial constriction of the vessel travelling along its length, similar to the oesophagus or the ureter [3,21]. Evidence supporting this view exists in the literature [32], with 84 per cent of the contractions observed propagating along the entire length of the lymphangion [33]. Peristalsis in guinea pig mesenteric lymphatics has also been reported [44]. The peristaltic character of contraction is disputed by others (J. E. Moore 2011, personal communication), and videos captured with a high-speed camera by Davis et al. [23] show the entire wall contracting at roughly the same time. It should be noted that these studies were conducted in different species: bovine in Ohhashi et al. [32] and rats in Davis et al. [23].
Lymph propulsion in collecting vessels is affected by preload, afterload, transmural pressure and shear stresses in addition to nerve and humoral mediators. These are termed intrinsic mechanisms. The details of these mechanisms and, especially, the molecular regulation are not fully understood [46]. Extrinsic mechanisms such as skeletal muscle contraction, motion of adjacent organs and arterial pulsations also influence lymph flow. External forces such as massaging, a common treatment modality for lymphoedema, were once thought to affect lymph flow but evidence shows that their effect is on the filling of initial lymphatics rather than on the pumping of the larger collecting vessels [47].
It is not clear whether extrinsic mechanisms can have a significant or even dominant role in the pumping of lymph. Olszewski [48] and Olszewski & Engeset [49] consider extrinsic mechanisms to be secondary to the intrinsic pumping based on their results in humans. On the other hand, Negrini et al. [50] found that, in rabbits, pleural lymph flow is due to respiratory movements of up to 40 per cent of the total volume, although it cannot be excluded that, at other sites, the intrinsic mechanism is not the dominant one.
Passive flow owing to a positive pressure gradient may also occur in oedema, during which lymph formation is increased [47]. In fact, Gashev et al. [51] report that, for low levels of lymph formation, the intrinsic mechanism is dominant but, as these levels rise, the active lymph pump is inhibited and the vessels become conduits. However, according to the authors, this does not prove that the flow-related inhibition decreases total flow in vivo, since, at higher levels of lymph formation, passive flow may be dominant.
The lymphatic structure and function are very complicated and can only be briefly addressed in the context of the present review. The reader is referred to the reviews by Aukland & Reed [40], Aukland [47], Gashev & Zawieja [22,41], Reddy [52], Schmid-Schonbein [34,43] and Swartz [3], the textbook by Földi et al. [7] and the references therein for more information. Unfortunately, the information on the human lymphatic system is scarce. Most studies have been performed in different animal species with different experimental procedures at different locations, which adds to the confusion. Extrapolation of results to humans should be performed with care as there are differences in species as well as differences between members of the same species. It is not surprising therefore that there are still structures in lymphatic vessels found in humans (ampullae and diverticula) for which no explanation of their function exists [17].
Review. Modelling the lymphatic system K. N. Margaris and R. A. Black 605

MATHEMATICAL MODELLING
Mathematical modelling of the lymphatic system is in its infancy. The first model of the entire lymphatic system is attributed to Reddy [53]. It is a simplified one-dimensional model of the entire lymphatic network that is based on the Navier -Stokes equations of fluid mechanics. Although parts of the lymphatic anatomy had to be excluded for simplicity or because of a lack of available information, this was a pioneering work at that time and remains the only model of the entire lymphatic circulation. In contrast, the arterial tree has attracted more attention, with upwards of 23 models cited in the literature to date [54].
Recently, the lymphatic system has attracted more attention. Macdonald  The models mentioned above are either lumped or one-dimensional. There are even fewer models in higher dimensions. In fact, to the authors' knowledge, only one three-dimensional model of a single lymphangion, albeit neglecting the valves, has been published so far by Rahbar & Moore [60].
A few models of the initial lymphatics have been developed as well. The first attempt is arguably that of Reddy & Patel [42]. Mendoza & Schmid-Schonbein [61] developed a two-dimensional model of the primary lymphatic valves, which was later refined by Galie & Spilker [62]. Recently, Roose & Swartz [63] developed a multi-scale model of the initial lymphatics using homogenization theory in order to study the optimal structure of the lymphatic capillary network with respect to interstitial fluid drainage.
Models of the secondary valves are currently lacking. Macdonald [64] included a very simple two-dimensional model of a secondary lymphatic valve, although this was a very crude model in which the motion of the valve leaflet was imposed rather than solved through a fluid -structure interaction model.
The models that exist so far in the literature can be classified into two categories: lumped and continuum models. In the context of the present text, a continuum model refers to a model described by partial differential equations as opposed to lumped models that are described by ordinary differential or algebraic equations. In § §4.1 and 4.2, the abovementioned models are discussed and their limitations are elucidated.

Lumped models
One type of approach that can be used to model lymph flow is the lumped approach in which a distributed system is reduced to a discrete one. Obviously, these models are simple and the mathematical complications are kept to a minimum. They can yield useful insight into the behaviour of the system under investigation; however, they cannot yield a detailed description of the flow.
Drake et al. [59] used an electrical analogue model of lymphatic vessels. The vessel was modelled as a resistance and a voltage source to take into account the flow resistance and pumping action of the vessel, respectively. A diode was added to account for the valves. This very simple approach allows the inclusion of a large number of vessels, or regions of the lymphatic circulation, by lumping them into an equivalent resistance and voltage and connecting the resulting equivalent series in series or parallel. However, the authors used the Thevenin theorem in their analysis, which imposes the requirement that the pressure-flow rate relationship is linear. Depending on the species or lymphatic compartment under consideration, this relationship is, in general, highly nonlinear [32,[65][66][67]. Thus, the applicability of such models is questionable.
The same group undertook additional work using a model composed of a linear circuit component and a nonlinear model of a lymphatic pump composed of five lymphangions. Both passive and active flow were taken into account and, despite its simplicity, the model showed that, when outflow pressure is higher than inflow pressure, flow is mainly due to the active contractions, but the pressure gradient becomes positive when most of the flow is due to passive factors [65,68]. This is in line with the experimental findings by Gashev et al. [51].
Quick et al. [57] developed a lumped model by employing the time-varying elastance approach, which was originally developed for the ventricles of the heart [69]. This model was compared with a more complex lumped model, derived by Quick and co-workers [70,71] using the linearized Navier-Stokes equations of fluid flow. Quick et al. [57] claimed that the simplification of the their model to a simple algebraic equation, by neglecting inertia and the viscosity effect, is justifiable and the error introduced in mean flow is small. However, in the vicinity of the valves, the inertial terms are expected to be dominant, according to Reddy & Kesavan [72].
Venugopal et al. [58] refined the time-varying elastance model using a bilinear approximation for the pressure -volume relationship, which is certainly an improvement as the actual P-V relationship is nonlinear and has been approximated with an exponential function by Ohhashi et al. [32]. The hypothesis proposed is that lymphangions change behaviour depending on the magnitude of the transmural pressure: at higher transmural pressures, lymphangions become insensitive to it and maintain constant stroke volume.
Bertram et al. [56] modelled a chain of lymphangions and concentrated on the pumping behaviour and derived pump-characteristic curves. As with the rest of the lumped models discussed in this section, a peristaltic wave of contraction cannot be modelled without the addition of more lumped components. This model added a pressure-dependent valve resistance, and is the first modelling attempt to do so. The constitutive equation of the vessel wall included passive elasticity and active contraction terms. The inclusion of these terms allowed for wall stiffening at high internal pressures and loss of compliance at high external ones. The resulting nonlinear ordinary differential equations yielded a complicated behaviour, which depended on the parameters and the number of lymphangions in a chain.
An interesting result of this model was that coordinated pumping between lymphangions is inefficient, which contradicts the results derived from the time-varying elastance model by Venugopal et al. [71]. The contraction waves have been observed in vitro to propagate downstream or upstream with a finite velocity of 5-8 mm s 21 . In vivo contractions can be irregular and lymphangions have a refractory period, which is not included in the work by Bertram et al. [56]. Lymphangions do not contract synchronously and consequently the result by Bertram et al. [56] is likely to be more realistic, even with the aforementioned limitations, than the ones obtained from the time-varying elastance models. In the work by Venugopal et al. [71], an experimental verification was also provided, but the experiments were conducted in valveless lymphangions. Therefore, the results or conclusions from this study must be interpreted with caution, owing to the importance of the one-way valve system in lymphatic flow generation. It was also found that pumping is less effective as the lymphangion diameter decreases owing to an increase in external pressure, indicating that when oedema is present collecting vessels in the affected region will not be able to contribute to lymph clearance by active contraction.
While simple approaches such as the work by Quick et al. [57] and Venugopal et al. [71] discussed in this section are useful, and can provide useful insight into the function of the lymphatic system when used with caution, their limitations should not be overlooked. The elastance is sensitive to pressure changes, and the nonlinear pressure-volume behaviour of lymphangions limits the applicability of this concept to the study of the lymphatic system [67]. Moreover, there are differences in the behaviour of lymphangions and ventricles [58] and the time-varying elastance model has been criticized as insufficient to describe the ventricle dynamics as well, for which it was initially developed [73,74].
Reddy [53] was the first to develop a model of a large part of the lymphatic network. The model is derived from the one-dimensional Navier -Stokes equations [72,75], by assuming that inertial effects are negligible. The lymphatic vessels were assumed to be contracting, although a simplification was made by assuming that the vessels have a uniform radius over the length, i.e. they retain a cylindrical shape. Spatial averages of the pressure were used for each lymphangion, and by doing so the Navier -Stokes equations reduced to ordinary differential equations, and hence this model is effectively a lumped one.
After deriving the equations describing fluid flow and momentum transport of a single lymphangion, Reddy and co-workers [55,75] proceeded to assemble them into a model of the thoracic duct, or a network model of the whole system, although several simplifications were necessary. The initial lymphatics were considered to be lumped. The nodes were not considered owing to the fact that information on the microcirculation inside them was not available, although it is thought that they affect lymph flow owing to their contractility, high resistance and mass exchange processes that take place inside them. The model was completed via the incorporation of the following boundary conditions: pressure in the jugular vein, at arterial and venous capillaries, capillary protein concentration and external pressure exerted on each lymphangion. It should be noted that this is the only model that takes into account the refractory period of lymphangions. The results, according to these authors, were in agreement with experimental observations in sheep by Hall et al. [76] when extrapolated in humans. However, the validity of this approach is questioned by Olszewski & Engeset [49].
There are several aspects of lymphangion behaviour that have not been addressed by lumped models; in fact, they are not addressed by continuum models either. These aspects are the dependence of contraction on transmural pressure and shear stresses, the influence of humoral factors and nerve stimulation.
Lumped models can nevertheless yield a very complicated behaviour, as shown by Reddy [53] and Bertram et al. [56], but cannot resolve space-dependent flow characteristics, at least without the inclusion of more lumped components. This is where more advanced computational models can have an advantage.

Continuum models
Continuum models of the lymphatic system are scarce. As mentioned previously, Macdonald et al. [45] refined Reddy's mathematical model by including two terms in the transmural pressure equation: a tension term and a damping term. Moreover, spatial averages of the pressure were not assumed, thus retaining the continuum nature of the governing equations. Two notable differences are that the contraction of the vessel wall was modelled by a change in Young's modulus rather than an active stress added to the stress state of the material, and it was assumed that lymphangions do not have a refractory period, as in the work by Reddy et al. [75]. Data for parameters such as the Young modulus and damping were derived by experiments on bovine mesenteric lymphatics, although the properties measured had a very large variation owing to uncertainties in vessel diameter and wall thickness measurements. 4 The Young modulus was found to be 1.2 + 0.7 kPa; a quite large variation, possibly even larger than the physiological variations that one would expect to find between individual members of species. The damping term also varied considerably.
The models discussed so far, either lumped or continuum, have some common limitations. First of all, they assume one-dimensional laminar flow with a parabolic velocity profile (i.e. fully developed), whereas in reality, based on anatomical findings, the length to diameter ratio varies widely depending on the anatomical site. While the laminar flow assumption may be valid for lengthy vessels, this is not necessarily the case for shorter vessels with valves at either end, where the flow field is likely to be substantially different. Reynolds numbers of up to 5 have been reported in the literature [77], and while these values are low it must be emphasized that, if the values reported are correct, the inertia forces are by definition five times greater than viscous ones, giving rise to the possibility of localized eddy formation [44] and areas of recirculation and mixing. (Even in the case of a Reynolds number of 1, the inertia and viscous forces are of equal magnitude and their effect should be studied further before concluding that they can be neglected.) Further complications may arise when the flow is unsteady, and the classical definition of Reynolds number, Re ¼ rud/m, which has been used in the literature to characterize lymph flow is not appropriate in this case. In peristaltic flows, the modified Reynolds number quoted by Shapiro et al. [78] is more appropriate: Re ¼ ra 2 c/ml, where c is the wave propagation velocity, a is the mean half-width of the channel (mean radius in cylindrical channels) and l is the wavelength.
While the formation of turbulence in the vicinity of the valves is possible according to Gashev [79], evidence for the existence of turbulence in lymphatics is not found elsewhere in the literature. Studies of flow within microfluidic devices show that transition to turbulent flow depends on the entrance length, and may take place at Reynolds numbers far lower than that of fully developed pipe flows [80]. Nevertheless, true turbulence is a phenomenon that needs time and space to develop, and the evidence available to us at this time suggests that its formation within the lymphatic system is highly unlikely.
Recently, Rahbar & Moore [60] challenged the assumption of one-dimensional flow by developing a three-dimensional computational model of a contracting lymphangion, without valves and inertial effects, and still assuming a slender vessel with a maximum aspect ratio of 3.75. This is not a true fluid -structure interaction simulation, though, as the wall motion was prescribed as a boundary condition.
Comparison of the wall shear stresses with those calculated from the Poiseuille equation yielded a difference of no more than 4 per cent in all the simulated conditions, despite the fact that the radial velocity can be 30 per cent larger than the axial velocity in certain cases. Poiseuille's law yielded a resonable estimation of the shear stresses; however, it remains to be validated in the case of a lymphangion with valves.
The treatment of the valves in these models is rather simplistic. Reddy et al. [55] assumed the resistance of the valves to be dependent only on the vessel diameter and imposed the condition that the flow rate is positive in all lymphangions. Essentially, these assumptions neglect the action (and possible interaction) of the valves. Macdonald et al. [45] assigned a timing of 20 ms, but with a constant resistance to flow for the one-dimensional model. Macdonald [64] also developed a computational fluid dynamics model of a valve, but the leaflet motion was imposed rather than calculated from the interaction with the flow. Bertram et al. [56] considered the valves to have a resistance which was a function of the pressure drop across the valve. However, the nature of this relationship in humans is unknown and data for other species are scarce as well.
Damage or removal of nodes increases the occurrence of lymphoedema when compared with vessel removal, thus indicating the importance of nodes in lymph flow [7]. Yet, none of these models includes the effect of lymph nodes on the propulsion of lymph. In fact, to the authors' knowledge, no modelling of the flow or mass exchange in the nodes has been attempted to date. Lymph nodes have smooth muscle in their walls and have been shown to contract with lower frequency than vessels. They are thought to act as reservoirs as well as pumps and filters. Their resistance to flow is quite high. Yet, there are a few studies of lymph node biomechanical behaviour; none of these includes lymph flow, but rather they concentrate on the contraction of the node's smooth muscle layer [29,81].
The computational models discussed to this point model the lymphatic system as a network of lymphangions; they are not concerned with the initial lymphatics. The first attempt to model the primary lymphatic valves should probably be credited to Reddy & Patel [42]. The flow was assumed to be axisymmetric and driven by a pressure difference between the interstitial space and the lymphatic vessel. The capillary was assumed to be supported by a linear elastic spring to account for the anchoring filaments. The authors added a hoop stress in the momentum equation of the vessel wall, even though the capacity of endothelial cells to produce a circumferential stress is not generally accepted [34]. Nevertheless, some interesting results were derived. It was found that, by increasing the amplitude of pressure pulses in the adjacent collecting vessel, the flow rate in the initial vessel increased; a result which is in line with the hypothesis that suction from collecting vessels facilitates lymph formation. Changes in interstitial fluid pressure also affected the filling of capillaries, which again is in agreement with the general concept of lymph formation. In general, though, the model suffered from lack of data regarding the initial lymphatics and therefore some parameters had to be estimated based on properties of blood vessels.
The optimal structure of the lymphatic capillary network was the focus of a study by Roose & Swartz [63]. Using homogenization theory, the authors developed a model that captures the influence of the microstructure of the capillary network on the macroscale fluid drainage at the tissue level. It was found that a hexagonal network is optimal in terms of fluid drainage, which is actually the structure found in mouse tail and human skin. It is interesting to note that the results obtained for the fluid clearance with respect to time at various positions along the network are qualitatively quite similar to those obtained by Reddy & Patel [42].
Mendoza & Schmid-Schonbein [61] developed a model of the initial lymphatics, but focused more on the endothelial cells that form the primary valves. The cells were modelled as elastic beams with the assumption of small deflections that were subjected to a fluid pressure difference generated by one-dimensional creeping flow into the lymphatics. It is interesting to note that, although 8 years had elapsed since the work by Reddy & Patel [42], the authors had to use data from vascular endothelial cells owing to the lack of data for lymphatics.
In the work by Mendoza & Schmid-Schonbein [61], small strains were assumed. The more elaborate model by Galie & Spilker [62] assumes large strains. This model also assumed two-dimensional flow in contrast to the model of Mendoza & Schmid-Schonbein [61]. The interstitial space was included as a porous medium, neglecting deformations and describing flow with Brinkman's equations. The model was completed by appropriate no-slip boundary conditions at the stationary walls, velocity matching at the fluid-structure interface, pressure and velocity continuity at the porous fluid interface and pressure at the boundaries. The authors concluded that these additions had an effect on the results as the calculated cell deflections were large. Both models still lack appropriate material properties for lymphatic cells.

FUTURE WORK AND CONCLUSIONS
The lymphatic system is a highly complex active system; a fact that when combined with the scarcity of anatomical and physiological data results in a formidable problem to model with mathematical methods, either analytically or numerically. In this section, a lumped model initially developed for the ventricle of the heart by Formaggia et al. [82] is proposed as an alternative model for a lymphangion. A lumped approach has the benefit of simplifying the problem, but at the same time yields a model that can be used to model large parts of or even the entire lymphatic circulation. The electrical analogue model proposed by Formaggia et al. [82] is shown in figure 5. Despite the fact that this was proposed for the heart, the resemblance of lymphangion function to that of a ventricle makes this model suitable as a starting point.
The lack of anatomical data favours the use of lumped models, as the spatial dimensions are reduced to zero. The lymphatic one-way valves are represented by diodes. The compliance of the wall is represented by a capacitor C and a time-varying resistance dC/dt. The action of the smooth muscle is modelled as a current source M(P,Q,t), which is a function of pressure P, flow rate Q and time t. The proposed model alleviates the limitations of the time-varying elastance concept used in some of the models discussed previously. It is also readily extensible. The addition of an extra resistor R 0 connected to the ground (shown as the dashed line box in figure 5) allows for leakage of the vessel to be taken into account. The addition of an inductor in series with resistor R (not shown) can account for inertial effects of the fluid.
The parameters of the components can be easily altered for parametric studies, or can be made a function of other factors such as neural or hormonal stimulation. Unit models can be connected to form chains of lymphangions and bifurcations in order to model larger parts, or even the entire lymphatic circulation. The model can be connected with lumped models of the systemic circulation as well as in order to study the interaction of the two systems. More complicated multi-scale models can be constructed, in which parts of the system are modelled as lumped models and others as continuum one-dimensional (or higher) ones, as has been done for the cardiovascular system [83]. Interestingly, the same model, with appropriate component parameters, can be used to model the lymph nodes as well.
The paradigm for modelling of physiological flows has been set by the effort undertaken in the cardiovascular system. The abovementioned model is not the only one that can be drawn from the heavily researched systemic circulation. The methodologies developed for the systemic circulation range from simple Windkessel models [84] and multi-physics/multi-scale [85] to models of cardiac cells [86] and excitable tissues [87], and in principle may be applied to the lymphatic system as well.
The models of the lymphatic system discussed in this paper as well as the models of the cardiovascular system discussed in the previous paragraph are all based on continuum mechanics, but this is not the only possible approach to modelling biological systems; in fact, a continuum-based approach is not the most suitable for describing mechanobiological events [88]. Alternatives to phenomenological continuum mechanics models that can incorporate the activities of cells also exist. Cellular automaton and agent-based models of cellular growth and cell -cell interaction have existed for quite some time in the literature. These methods have also been coupled with continuum mechanics to study arterial adaptation [88,89]. Another approach is the theory of active particles [90][91][92][93]. It extends the kinetic theory of particles by assigning a state that defines the microscopic behaviour of each particle. Thus, macroscopic models of tissues can be derived based on the underlying behaviour of cells, cell -cell interaction, cell proliferation and death. These powerful mathematical tools afford the means to integrate events at the molecular, cellular and tissue levels, with the prospect of improved models which will find application in developmental biology, angiogenesis, lymphangiogenesis, tissue remodelling and tumour growth as well as interaction of cancer cells with those of the immune system.
In relation to the cardiovascular system, in which modelling efforts have been numerous and have advanced to the level of multi-scale and multi-physics even including models of excitable tissues, the work carried out on the lymphatic system seems to be in its infancy. More work will be necessary in combination with experimental verification in order to progress the knowledge on the function of the lymphatic system in health and disease, but, as our knowledge and Review. Modelling the lymphatic system K. N. Margaris and R. A. Black 609 understanding of the lymphatic system increases, so too will the likelihood that new treatments will emerge.