A Mathematical Description of the Flow in a Spherical Lymph Node

The motion of the lymph has a very important role in the immune system, and it is influenced by the porosity of the lymph nodes: more than 90% takes the peripheral path without entering the lymphoid compartment. In this paper, we construct a mathematical model of a lymph node assumed to have a spherical geometry, where the subcapsular sinus is a thin spherical shell near the external wall of the lymph node and the core is a porous material describing the lymphoid compartment. For the mathematical formulation, we assume incompressibility and we use Stokes together with Darcy–Brinkman equation for the flow of the lymph. Thanks to the hypothesis of axisymmetric flow with respect to the azimuthal angle and the use of the stream function approach, we find an explicit solution for the fully developed pulsatile flow in terms of Gegenbauer polynomials. A selected set of plots is provided to show the trend of motion in the case of physiological parameters. Then, a finite element simulation is performed and it is compared with the explicit solution.


Introduction
Lymph nodes are organs scattered throughout the lymphatic system which play a vital role in our immune response in breaking down bacteria, viruses and waste; the B Alessandro Musesti alessandro.musesti@unicatt.it Giulia Giantesio giulia.giantesio@unicatt.it Alberto Girelli a.girelli3@campus.unimib.it 1 Dipartimento di Matematica e Fisica "N. Tartaglia", Università Cattolica del Sacro Cuore, Brescia, Italy 2 Dipartimento di Matematica e Applicazioni, Università degli Studi di Milano-Bicocca, Milan, Italy interstitial fluid (called lymph once inside the lymphatic system) is of fundamental importance in doing this since it transports these substances inside the lymph node (Arasa et al. 2021). The main features of the lymph node from a mechanical point of view are the presence of a porous bulk region (lymphoid compartment, LC), surrounded by a thin layer (subcapsular sinus, SCS) where the fluid can flow freely. More than 90% of the lymph remains in the SCS, while the remaining part enters into the LC through a conduit system network (Roozendaal et al. 2008;Grebennikov et al. 2016;Savinkov et al. 2017) formed by fibroblastic reticular cells (FRC), which is the parenchyma of the LN (Novkovic et al. 2020); due to this, LNs are organs with high resistance to flow. The bigger particles cannot enter the conduits formed by FRC and remain in the SCS, where they are confined and filtered by specialized cells of the lymphatic endothelial cells; however, there is some evidence that the selectivity of the FRC network is not based solely on the size of the molecules; indeed, selected macromolecules, such as antibodies, can gain access to the LN parenchyma (von Andrian and Mempel 2003).
Lymph flow inside LNs has an important function; indeed, fluid flow biases macromolecular distribution, enhances ligand expression, aligns extracellular matrix and shapes active mechanisms of cell migration. Fluid flow through endothelial monolayers and FRC networks enhances the expression of chemokines that direct leucocyte localization and migration patterns (O'Melia et al. 2019). Moreover, increased flows enhance proliferation and drug sensitivity in B cell lymphoma (Apoorva et al. 2018;Lamaison et al. 2020). Fluid flow is important to study the tumor metastasis (Birmingham et al. 2020) and drug transport (Permana et al. 2021). Despite its importance, as far as we know, only few models in the literature try to describe the behavior of the lymph from a mechanical point of view (Novkovic et al. 2018;Jafarnejad et al. 2015;Cooper et al. 2016Cooper et al. , 2018Tretiakova et al. 2021;Giantesio et al. 2021) or mimicking the LN mechanical properties in a LN-on-a-chip model (Shanti et al. 2020;Birmingham et al. 2020;Shanti et al. 2018).
In this paper, we propose a mathematical model for the flow of the interstitial fluid in a lymph node. We assume the lymph to be an incompressible fluid similar to water; moreover, we assume a small Reynolds number as a result of the small velocities within the lymph nodes (Moore and Bertram 2018), hence we can model the flow into the LC by Darcy-Brinkman equation [due to the high porosity and the time-dependence of the flow (Shanti et al. 2020;Savinkov et al. 2017)], and the flow inside the SCS by Stokes equation. The lymph enters the lymph node from the lymphatic vessels, which have a complex structure formed by one-way valves that prevent retrograde flow and a wall structure composed of sinus-lining cells: such cells control and generate active pulsation of the wall, pumping the lymph from a segment between two valves to another (the segment is called lymphangion) (Mozokhina and Savinkov 2020;Moore and Bertram 2018). This means that the lymph has a relevant pulsatile behavior, and we take it into account in our model. In Sect. 2 we describe the behavior of the lymph explicitly in spherical geometry, supposing that the fluid flow inside the lymph node is axisymmetric with respect to the azimuthal angle, so that we can assume a simplified two-dimensional geometry and we can use the stream function approach (Happel and Brenner 1983) to find an explicit solution. We remark that the solution given in Sect. 2.1 is quite general and can be used also for other choices of boundary conditions. Finally, in Sect. 3 we compare our results with some finite element simulations obtained using the open source software FreeFEM (Hecht 2012).

Explicit Result in a Simplified Case
Let us model the lymph node (LN) as a spherical region: the subcapsular sinus (SCS) is a thin spherical shell with radii R 1 < R 2 of creeping fluid flowing near the external wall of the LN, while the lymphoid compartment (LC) is a sphere of radius R 1 of porous material. We use spherical coordinates (r , θ, φ), where r is the radial distance, θ the polar angle and φ the azimuthal angle; moreover, we suppose axial symmetry with respect to the azimuthal angle φ.
Assuming that the lymph, which flows inside the LN, is an incompressible fluid, and that the Reynolds number is small, we have the equations where ρ 0 is the constant density, v the velocity, p is the pressure, μ the viscosity of the lymph, μ e the effective viscosity, k the permeability. The second equation in (1) is the Stokes equation and describes the motion in the subcapsular sinus, the first is the Darcy-Brinkman equation, which is used for modeling the flow in the porous region of the LC, while the last equation models the incompressibility of the fluid.
Here we assume a constant homogeneous permeability k (Savinkov et al. 2017;Shanti et al. 2020). The effective viscosity μ e in general differs from the classical viscosity μ because μ e keeps into account the Brinkman correction (Nield 2000). Furthermore, assuming that the flow is time periodic with period T , we write the time dependence of the velocity and of the pressure as a Fourier expansion where ω = 2π/T .

Solving the Equations
Now we want to compute the general solution of system (1) in terms of the Fourier expansion (2). Here we try to be as general as possible, without imposing any boundary condition, so that out solution can be used in several situations. We will deal with suitable boundary conditions for our specific problem in Sect. 2.3.
By using (2), system (1) becomes which can be written in compact form as where Z is the set of integers, while q m is given by Moreover, it is useful to perform the change of variable ζ := cos θ , so that the previous equations become By introducing the operator we can rewrite (4) as while for the pressure we have Focusing on the case m = 0, we have that the solution can be written as We can now solve (10): by using the separation of variables substituting in the first equation of (10) we get As the first term of (12) depends only on r and the second term only on ζ , the two have to be constant, say n(n − 1) with n ∈ N (Haberman and Sayre 1958), where N is the set of natural numbers. Hence (12) becomes The solution of (13) is given by for some constants A (n) , B (n) , while (14) is the Gegenbauer equation, whose solutions are the Gegenbauer polynomials G n , H n with order −1/2, of the first and second kind, respectively. Hence the solution of the first equation of (10) becomes m . Since H n is not smooth in ζ = ±1 and G 0 , G 1 lead to an infinite tangential velocity, the solution simplifies as m . The second equation of (10) is and, using again the separation of variables, by a similar procedure as before, we obtain Equation (18) is a Bessel equation, hence the solution can be written as where J s , Y s are the Bessel functions of the first and second kind, respectively. Equation (19) is the same Gegenbauer equation as (14), hence the solution of (17) is given by Now we want to employ the definition of q m , so that we have to distinguish between the Stokes and the Darcy-Brinkman case. Let us denote with A m those of the Darcy-Brinkman case (0 ≤ r ≤ R 1 ). Using (5), we obtain, for any m = 0, where the superscript S denotes the Stokes case and B the Darcy-Brinkman case, and we used the fact that r = 0 is in the domain of ψ B , so thatB Regarding the pressure, we use (9) to obtain in the Stokes case, and in the Darcy-Brinkman case, where P n are the Legendre polynomials of the first kind. For m = 0 we get the well-known steady solution of the Stokes equation and for the Darcy-Brinkman equation we have

Geometrical and Physiological Parameters
We use an idealized spherical geometry based on the data obtained from a murine (popliteal) lymph node: the radius is R 2 = 0.5 mm, the subcapsular sinus (SCS) thickness is h = 10 μm, the afferent and efferent lymphatic vessels have the same radius R LV = 40 μm (Birmingham et al. 2020;Kislitsyn et al. 2015;Ulvmar et al. 2014;Das et al. 2013;Shanti et al. 2020;Zhang et al. 2013;Jafarnejad et al. 2015).
With these data, we have that more than 90% of the lymph takes the peripheral path without entering the LC in a pulsation cycle (Jafarnejad et al. 2015;Guyton 1983, 1985). The inlet and outlet conditions are imposed in the upper and lower lymphatic vessel (near θ = 0 and θ = π , respectively) as a pulsatile flow of the form where L is the maximum lymph mean flow of the inlet lymphatic vessel. Here we assume L = 10 −3 mm 3 /s, as measured in (Blatter et al. 2016), and f (t) is a periodic function. The function H is given by where the constant 0 < ζ 0 < 1 describes the inlet and outlet regions, and is given by Notice that we are assuming that the inlet and outlet velocities are the same. The lymph is modeled as an incompressible Newtonian fluid similar to water (Moore and Bertram 2018) with viscosity μ = 1 mg/(mm s) and density ρ 0 = 1 mg/mm 3 . The permeability is considered homogeneous (Savinkov et al. 2017) with value k = 3.84 × 10 −9 mm 2 (Shanti et al. 2020). The effective viscosity is taken as μ e = μ φ (Ochoa-Tapia and Whitaker 1995a; Tan and Pillai 2009), where φ is the porosity taken as φ = 0.75 (Shanti et al. 2020). The parameters are summarized in Table 1.

Boundary Conditions
We now want to impose suitable boundary conditions to our general solution. We give a Dirichlet boundary condition at the external boundary and the Ochoa-Tapia boundary conditions (Ochoa-Tapia and Whitaker 1995a, b) at the interface between the porous zone LC and the free-fluid region SCS. In this way we can close the problem and find a unique solution. More precisely, we will assume the no-slip condition for the velocity on R 2 , except near θ = 0, π, where we impose the inlet/outlet flow (27). For simplicity, given the small diameter of the afferent/efferent lymphatic vessel, we impose the inlet/outlet condition only for the radial velocity v r , but we could use the same procedure to impose boundary condition for v θ too. For the boundary conditions on the internal radius R 1 , the Ochoa-Tapia boundary conditions imply the continuity of radial and tangential velocity, the continuity of the normal stress tensor and a jump-condition on the shear stress.
Thanks to the above conditions, we can determine for every n the six unknown constants in Eqs. Expanding the step function H (ζ ) in (28) in terms of Legendre polynomials, we get and we kept into account that b 0 = 0 since H is an odd function.
To impose the boundary condition, we need to expand in Fourier series the time dependence of (27), as we did in (2). Writing

Now we impose the boundary condition
for any m ∈ Z and n ≥ 2, where we used the linear independence of the Gegenbauer polynomials. By the no-slip boundary condition on v θ , recalling (7) 2 it follows that where we used again the linear independence of the Gegenbauer polynomials. We now write in terms of the stream function the Ochoa-Tapia boundary conditions on the internal radius R 1 (Prakash 2020), using the linear independence of Legendre and Gegenbauer polynomials: • Continuity of v θ : • Continuity of normal stress: where T rr,m = −p m + 2μ ∂v r ,m ∂r ; we can write this condition as • The stress jump condition: where β is the slip constant which has to be estimated experimentally. Since the expression of the shear stress is in the term ∂v r ,m ∂θ = − 1 − ζ 2 ∂v r ,m ∂ζ there is a second derivative of the Gegenbauer polynomials, so that we need the following property (Abramowitz and Stegun 1964;Zlatanovski 1999): Hence we have: After some computations, Eq. (37) can be written as From (31) Moreover, we fix the value of the pressure in one point to find the constants in equation (23)-(25) and have a physiological pressure value. By (36), it follows that C S m = C B m and C B 0 = C S 0 . We fix the pressure (with respect to time) at the exit point (r , ζ ) = (R 2 , −1) by using the same time function of (27), that is, Hence we can find the pressure constants by imposing where p S m (r , ζ ) is given in (23) for m = 0, and in (25) 2 for m = 0.

Explicit Results
This section is devoted to show some plots related to the explicit solution and to make some considerations about the proposed model. Following (Bertram et al. 2017), we choose a time function of the form We notice that in this case the period of a pulsatile flow in the lymph node is 2 s, hence ω = π , and f m = 0 for m = −1, 0, 1. In this model we do not take into account the inhibition and the autoregulation of the contractions in the lymphangion, given by several factors like shear stress and pressure (Bertram et al. 2019; Moore and Bertram 2018); a further extension of this model can be the coupling with a lymphangion model for taking into account these phenomena.
In Fig. 1, we plot the pressure distribution in the LN with the fixed constantp = 6.18 × 10 5 mPa, corresponding to the lower limit of the pressure found in Bouta et al. (2014); as we can see, the values of the pressure belong to the range given in that paper and, due to the incompressibility of the flow, the pressure translates from a higher value in the inlet zone to a lower value in the outlet zone. We can choose to fix any pressure at the outlet, and we have the same pressure distribution with different values, for example with the fixed pressure ofp = 4 × 10 5 mPa ≈ 3 mmHg as in Jafarnejad et al. (2015). Figure 2 provides the Stokes shear stress given by the formula: (in mPa) at time t = 1 s, where we have the maximum value of the velocity (and, consequently, of the shear stress) and radius r = R 1 (this is the shear stress at the exterior of the LC). We plot the shear stress value with two different boundary velocities: v in ≈ 0.22 corresponds to the physiological value of L = 10 −3 mm 3 /s, given in Table 1, found in Blatter et al. (2016), and v in ≈ 0.58 appears in Jafarnejad et al. (2015). As we can see, the shear stress is similar to the one reported in Birmingham et al. (2020) and Jafarnejad et al. (2015); that is, higher near the inlet flow and lower near θ = π 2 . The same behavior occurs in the velocity too (see Fig. 3). This trend is interesting because the cell adhesion to the exterior of the LC is proportional to the shear stress (Birmingham et al. 2020), hence the majority of the cells adhere (and then enter in the LC) near the inlet zone of the lymphatic vessel. Indeed, in our model the inlet shear stress is the same as the outlet one due to the choice of the same inlet/outlet velocity and the incompressibility of the fluid; however, usually a part of the lymph  (r , θ, t) in mPa with respect to the polar angle (θ = 0 near the inlet flow and θ = π near the outlet flow) calculated at t = 1 s and in the internal radius R 1 with different boundary velocities in mm/s (where v in ≈ 0.22 corresponds to L = 10 −3 mm 3 /s and v in ≈ 0.58 corresponds to L = 2.2 × 10 −3 mm 3 /s) (color figure online) enters in the blood capillaries in the LC Guyton 1983, 1985), so that the shear stress in the outer zone reduces.
As we can see in Figs. 3 and 4, for θ > 0 the tangential component v θ of the velocity in the SCS is the larger one. From the first picture in Fig. 3 one can see that the fluid flow in the porous medium is flat and starts increasing near the interface that connects the LC to the SCS, showing a non-differentiable point due to the Ochoa-Tapia boundary conditions (indeed, we do not impose the continuity of the derivative of v θ ).

Numerical Simulation
The explicit model found in the previous section uses several simplifications. In this section we propose some numerical simulations to describe a more general fluid flow in a lymph node.
We  The weak formulation of our problem is (supposing a constant density ρ = ρ 0 and viscosity ν = μ/ρ 0 ): for all w ∈ W S g , w b ∈ W B g such that w = w b on , and for all q ∈ Q S and q b ∈ Q B . In equation (42) we have that v is the velocity in S , p ∈ Q is the pressure in S , v b is the velocity in B , p b ∈ Q is the pressure in B , D(v) = 1/2(∇v + ∇v T ), T = −pI + μ ∇v + ∇v T , ν e = μ e /ρ 0 , K is the permeability tensor (in the case of Sect. 2, K = kI), T e = −p b I + μ e ∇v b + ∇v T b . Now we want to write the weak formulation for the boundary condition 2.3; we have that the continuity of the velocity is verified automatically, and, for the stress-jump condition, we have (Tan and Pillai 2009) on the interface : where B is the slip tensor (in the case of Sect. 2.3, B = βI).
The boundary conditions in the external wall (inlet condition and no-slip boundary condition) are imposed by the penalty method. Moreover, we add the Grad-div stabilization terms in either Stokes and Darcy-Brinkman domains (Jenkins et al. 2014;Neilan and Zytoon 2020;Qin et al. 2020;Rong and Fiordilino 2020). Thanks to this stabilization term, we have the stability for the Darcy-Brinkman equation (see Xie et al. 2008). For the numerical discretization, we use a BDF2 method for the time discretization, instead, we use P d k −P k element pairs (where k is the polynomials order and d is the dimension) with the Brezzi-Pitkäranta stabilization, which consists in adding the term S ∇ p · ∇qdV + B ∇ p b · ∇q b dV to the discretization of the Eq. (42), with ≈ h 2 T , where h T is the maximum diameter of the triangle of the finite element triangulation. The weak formulation here proposed has been implemented using the open source software FreeFEM.

Numerical Test
In this section we want to qualitatively compare the results obtained with the numerical simulation with the explicit results exposed in Sect. 2. For that, we use the same geometry and parameters exhibited in Sect. 2.2; hence, in the external boundary we will impose only Dirichlet boundary condition ( N is empty), subdivided as D = in ∪ out ∪ BC , where we are imposing the inlet and the outlet flow in in and out , respectively, given by the Eq. (27) with L = 10 −3 mm 3 /s, and the no-slip boundary condition in BC . The numerical stabilization parameters are estimated as γ 2 = 0, while γ 1 = 300 in S and γ 1 = 10 6 in B .
In Figs. 5, 6 and 7, we can see the tangential and radial velocity, and the shear stress, respectively. We can see that the results are very similar to the ones explicitly in Sect. 2.4: in order to remove some small oscillations in the internal velocity near R 1 , we needed to use a finer mesh, which meant a greater computational cost for every time step. We can do only a qualitative comparison between the numerical solution of this section and the explicit solution in Sect. 2.4 because we do not have available and precise physiological data of the lymph node and we have an error in both cases: in the explicit result from the truncation of the sum, and here due to the finite element approximation. Qualitatively, we have the same behavior and values here and in the explicit result.

Numerical Results
In this section we want to show a more complete numerical simulation using the method given in Sect. 3.
We use a spherical idealized 2D geometry with the same parameters given in Sect. 2.2; hence we suppose that the permeability tensor K is homogeneous and con- Fig. 5 Tangential component of the velocity in mm/s with respect to the radius at different angles at t = 1 s. The first graph corresponds to the tangential velocity in the LC (porous part), and the second corresponds to the tangential velocity in the SCS (free-fluid region) (color figure online) stant (this is not a limiting assumption, see Savinkov et al. (2017), Shanti et al. (2020), andJafarnejad et al. (2015)) and the same with the slip tensor B = βI. Moreover, we add to the simulation domain a part to the inlet and outlet lymphatic vessel (see Figs. 10 and 8).
As we mention in Sect. 2.2, more than 90% of the lymph takes a peripheral path; the lymph that enters in the LC does not remain in the LC but gets out due to the incompressibility of the lymph, because we are not taking into account the fluid exchange behavior given by the blood vessels inside the LN.
The inlet condition is imposed in the upper lymphatic vessel as a uniform pulsatile flow in the y direction with the Eq. (27). For the outlet condition, we need to fix the stress. For clarity and for a simpler interpretation, we fix the pressurep(t) in this way: We use the numerical parameters given in Sect. 3.1. In Fig. 8, we can see the pressure distribution in the LN withp = 6.18 × 10 5 f (t) mPa (= 6.3 cmH 2 O as the inferior limit in the range of pressure found in  Bouta et al. (2014) and as in the explicit results in Sect. 2.4), where f (t) is the one given by the Eq. (41). As we can see, the pressure distribution is similar to the one in Fig. 1 and it is in range with the corresponding results. If one hasp = 4×10 5 f (t) mPa (= 3 mmHg as in Jafarnejad et al. (2015)), the behavior of the pressure is similar to the one showed in Fig. 8 (so that we omit the picture), with a range of values comparable to Jafarnejad et al. (2015).
In Fig. 9, we can see the shear stress over time (in mPa). At time t = 1 s, we have the maximum value of the velocity (and, consequently, of the shear stress) and the shear stress is similar to the one found in the explicit result (the blue curve with v in ≈ 0.22 plotted in Fig. 2), that is in range with the values found in Birmingham et al. (2020) and Jafarnejad et al. (2015). We can see the norm and the velocity behavior in more details in Fig. 10. The tangential velocity (the most relevant one) is shown in Fig. 11. As expected, the maximum velocity is in the SCS near the inlet and the outlet region. In particular, we can see that the maximum velocity is between the connections of the SCS with the afferent/efferent vessel; then the velocity decrease with respect to the polar coordinate θ , reaching the minimum at θ = π/2. Moreover, even if we do not impose the outlet velocity equal to the inlet one, we have that this is true due to the incompressibility; hence our assumption used to find the explicit solution is not too limiting in this case.

Fig. 11
Tangential component of the velocity in mm/s with respect to the radius at different angles at t = 1 s. The first graph corresponds to the tangential velocity in the LC (porous part), and the second corresponds to the tangential velocity in the SCS (free-fluid region) (color figure online)

Conclusion
We proposed a model that describes the pulsatile lymph flow inside a simplified spherical lymph node (LN), using the Darcy-Brinkman equation to describe the lymph flow in the lymphoid compartment (LC, the porous part) and the Stokes equation to describe the flow inside the subcapsular sinus (SCS, the free fluid region). We found the explicit solution in terms of Gegenbauer polynomials and we showed the trend of the velocity, the pressure and the shear stress inside the LN; after that, we used this explicit solution to validate the numerical simulations of the model. Finally, we performed a more general numerical simulation with finite elements.
This model allows to better understand the fluid behavior inside the LN and how it changes with respect to time. The results obtained by our model are in agreement with the literature (Birmingham et al. 2020;Jafarnejad et al. 2015;Shanti et al. 2020;Cooper et al. 2016Cooper et al. , 2018. We remark that the Ochoa-Tapia boundary condition minimally affects the fluid behavior in the SCS. Still, it affects the flow in the LC, inducing a velocity profile which is not smooth at the interface between the LC and SCS regions. Particular attention was paid to the shear stress, because a lot of biological phenomena in the LN depend on it. Among them is the cell adhesion to the exterior of the LC, which is proportional to the shear stress: this is important because inside the LC there is a connection between the lymphatic system and the blood system, and some cells can get access from here to the blood circulation [for instance, tumor cells (Birmingham et al. 2020)]. Moreover, shear stress drives drug delivery and can affect pathologies like B-cell lymphoma (Apoorva et al. 2018;Lamaison et al. 2020). In our model we found that the shear stress is higher near the inlet and the outlet regions, and decreases with respect to the polar angle θ , reaching the minimum at θ = π/2; hence we believe that the majority of the cell adhesion is located near these two critical regions, which are the connections of the SCS with the afferent/efferent vessels.
Let us now make some considerations that can be interesting to improve the model in future. Here we have proposed to use a spherical geometry, but in general the LNs have a spheroidal shape (Giantesio et al. 2021;Jafarnejad et al. 2015;Cooper et al. 2016Cooper et al. , 2018. For simplicity, we did not take into account the fluid exchange inside the LN between the lymph in the fibroblastic reticular cells FRC and the blood in the capillaries, although it is important for the fluid regulation of the LN (Tobbia et al. 2009); a further extension of this work could take this phenomenon into account.
Moreover, in order to close our model and to find the unknown constants in the explicit solution, we used the Ochoa-Tapia boundary conditions, although other boundary conditions can be taken into consideration. For instance, a common choice is to impose the continuity also of the shear stress at the interface, which is tantamount to choose β = 0 in Eq.(40). Using the same technique, one can address the more general conditions given in Angot et al. (2017) and Angot (2018), where there is a discontinuity also of the tangential velocity and the normal stress.
Another interesting and important extension of this model would be to couple the flow in the lymph node with the flow in the lymphangion, in order to have more realistic time pulsation and see how the lymph node affects and regulates the global lymph circulation (Bertram et al. 2017(Bertram et al. , 2019Moore and Bertram 2018).
in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.