1 Introduction

The immune response requires recruitment of leukocytes to defend the body against foreign microorganisms such as bacteria and viruses. The leukocytes first form weak adhesion and roll on the endothelium surface. Then they bind firmly on the endothelial surface, pass through the endothelium in a process called diapedesis and breach the basement membrane. Finally, they move toward the chemotactic stimulus in the tissue [1]. Thus, reaction to infections requires frequent crossing of leukocytes through the layer of endothelial cells (ECs). To initiate the transmigration, leukocytes extend invadosome-like protrusions (ILPs) into the endothelium [2]. They can transmigrate either directly through the body of an individual EC (the transcellular route) or through the junction between ECs (the paracellular route) [2]. It is unclear what determines the path of transmigration.

Fig. 1
figure 1

Schematic of a leukocyte extending invadosome-like protrusions (ILPs) on the endothelium into the EC cell body and through the EC cell junction. The latter consists of transmembrane adhesive proteins (e.g., VE-cadherins) bridging a gap between the two ECs

Full size image

To address this question, Martinelli et al. [3] conducted three in vitro tests to examine the correlation between junctional integrity and the route of diapedesis. First, they compared the diapedesis through rat brain and heart and human heart and lung endothelia. The rat brain endothelium has stronger junctional integrity than heart and lung endothelia. They observed that for the rat brain endothelium the number of transcellular diapedesis is higher than the paracellular diapedesis, while in rat heart and also human heart and lung endothelia, the usage of paracellular route is dominant. Second, they used drugs and hormones to enhance or disrupt endothelial junctions in rat brain and heart tissues. They noticed that disrupting junctional integrity leads to a remarkable increase (about two folds) in paracellular diapedesis accompanied by decrease in transcellular transmigration. Finally, they exposed human lung endothelium to long-term shear flow as a mechanical modifying agent to promote the junctional strength and remodeling of the cytoskeleton. This alteration causes a significant increase in transcellular diapedesis. Based on these observations, Martinelli et al. [3] concluded that strong junctional integrity is correlated with dominant transcellular route of transmigration, and the EC junction tightness and local stiffness are the major determinants of the route of diapedesis. They hypothesized that leukocytes choose a path with the least mechanical resistance during transmigration, a tendency termed tenertaxis [3, 4].

Conceptually, the hypothesis appears plausible. But from a physical viewpoint, it raises several interesting questions. For example, as the leukocyte extends ILPs into the EC, how much resistance can the EC generate in either the transcellular or paracellular route? Given the typical level of protrusion forces in the ILPs, will they be able to overcome the resistance to effect a transmigration through either route? Will the difference in resistance correspond with experimental observations of the prevalence of one route over the other, under wild-type conditions and active intervention?

In this study, we test the tenertaxis hypothesis quantitatively by using a computational model of the endothelium and leukocytes with ILPs. To answer the questions raised above, we perform two subtasks. First, we review existing data on the mechanical properties of ECs in different tissues in their physiologic or altered states, as well as data on the shape and size of the invadosomes and representative magnitudes of the protrusive forces inside them. Thus, we identify a reasonable range of the mechanical parameters corresponding with prior experimental observations. Second, we use continuum theories of hyperelasticity and contact mechanics to simulate the penetration of an ILP through an EC body (transcellular) and through an EC–EC junction (paracellular). Our model predicts lower resistance in one route or the other depending on the endothelial properties in its natural or altered state. These outcomes are consistent with experimental observations of the prevalence of trans- or paracellular diapedesis. Thus, the model tentatively confirms the tenertaxis hypothesis.

2 Problem setup and methodology

According to the tenertaxis hypothesis, leukocytes extend ILPs into the EC to seek the route of least resistance [4]. Essentially, this amounts to comparing the mechanical resistance of the two potential routes: one against the formation of a tunnel through the EC body and the other against a protrusion through the endothelial junctional opening. Therefore, our model is centered on computing the mechanical resistance of the ECs monolayer against the protrusive force of the leukocytes in either route. As such, the model needs to represent at least the following three components (Fig. 1):

Fig. 2
figure 2

Geometric setup for the simulations of a the paracellular protrusion and b the transcellular protrusion. Initially, the tip of the rod is 10 nm above the flat apical surface of the endothelium in (a) and 20 nm in (b)

Full size image

(a) The ILP. We model the leukocyte protrusion as a nearly rigid rod with a diameter of 340 nm based on the average of protrusions observed in experiments [2, 4]. To simulate the rod penetrating the EC cell body or the EC junction, we can either specify the force on the rod or prescribe its displacement. Most of the results will be presented according to the latter protocol for simplicity and convenience.

(b) The ECs and the basement membrane underneath. We treat the ECs and the basement membrane as hyperelastic solid components obeying the neo-Hookean constitutive equation, ignoring the cytosol and other organelles inside the cell [5]. The EC nucleus is very stiff and diapedesis rarely happens in the area around the nucleus [4]. Therefore, we do not consider penetrations through the nuclear bulge of Fig. 1.

(c) Endothelial cell–cell junctions. We model the EC junction as a preexisting narrow gap between two neighboring ECs, with an undeformed gap size of 10 nm [6, 7]. The gap is bridged by bonds that consist of transmembrane adhesive proteins such as vascular endothelial cadherins (VE-cadherins) [8]. We model the molecular bonds as linear springs distributed between two adjacent EC edges [9, 10].

2.1 Geometric setup

Figure 2 shows the two geometries used for simulating the penetration of an ILP through the EC junction (Fig. 2a) and through the body of an EC (Fig. 2b). Because the diameter of the ILP is only about 2% that of the EC [4, 11], the penetration is largely a local event that does not involve the entire EC body. Thus, we circumscribe a small portion or subdomain of the endothelium for the simulations to reduce the computational cost. The subdomain used for the modeling of the paracellular protrusion is a three-dimensional cylinder 3 \(\upmu \)m in diameter and includes two adjacent ECs, the junction and the basement membrane underneath (Fig. 2a). The undeformed heights of the EC and the basement membrane are, respectively, 1 \(\upmu \)m and 2 \(\upmu \)m [12, 13]. The junctional gap has an undeformed width of \(d_{gap}=10\) nm [6, 7]. Its top opens into a groove, with two opposed circular arcs of radius 200 nm subtending a central angle of \(90^\circ \). For the transcellular protrusion, we use a two-dimensional axisymmetric domain 1.5 \(\upmu \)m in radius that comprises a part of the EC body (undeformed height of 1 \(\upmu \)m, away from the nucleus and the junction) and the basement membrane underneath (undeformed height of 2 \(\upmu \)m) (Fig. 2b). The protrusion is modeled as a rod 340 nm in diameter [4], with a cylindrical body and a hemispherical head.

To simulate the penetration process, we prescribe the kinematics of the movement of the ILP into the EC in terms of its tip displacement \(d_\mathrm{rod}(t)\). As the EC deforms, the contact force on the ILP is computed by integrating the traction over the ILP surface in contact with the EC. This yields the resistance force on the ILP, \(F_{r}(d_\mathrm{rod})\), as a function of the depth \(d_\mathrm{rod}\). Because the elastic deformation happens instantaneously, without viscous damping, essentially we are computing a series of quasi-static states with the rod at different positions of insertion. Thus, the speed of motion \(d_\mathrm{rod}'(t)\) plays no role in the result. In addition, we have carried out dynamic simulations by specifying a constant pushing force on the ILP and tracking its movement in time. The two protocols essential confirm each other, and we will mostly present data on the quasi-static \(F_{r}(d_\mathrm{rod})\) for simplicity.

The simulation ends when the displacement \(d_\mathrm{rod}\) reaches a threshold \(d_b\) for breakthrough. This threshold is set to the thickness of the EC monolayer (\(d_b=\)1 \(\upmu \)m) for the paracellular protrusion and \(d_b=0.98\) \(\upmu \)m for the transcellular protrusion. The threshold of 0.98 \(\upmu \)m for the transcellular protrusion was chosen because the thickness of the EC membrane is 10 nm [14, 15], and at this point, the apical membrane reaches the basal membrane so they may open up a transcellular tunnel [4]. After breakthrough, the leukocyte tends to spread between the EC and the basement membrane in vivo [16] or a substrate in vitro [3, 4], and we will not model that spread.

2.2 Governing equations and numerics

We treat the ECs and the basement membrane as neo-Hookean hyperelastic materials that can capture the strain-stiffening behavior of biological materials [17], both obeying the following constitutive equation with different material properties:

$$\begin{aligned} \varvec{{\varvec{\sigma }}}= G \, J^{-\frac{5}{3}} \left( \varvec{{\varvec{F}} {{\varvec{F}}}}^{T} - \frac{{\varvec{I}}}{3} \hbox {tr}(\varvec{{\varvec{F}} {{\varvec{F}}}}^{T}) \right) - K (J-1) \varvec{{{\varvec{I}}}}, \end{aligned}$$
(1)

where \(F_{ij}=\partial x_i/\partial X_j\) is the deformation gradient tensor, with \(\varvec{{\varvec{X}}}\) and \(\varvec{{\varvec{x}}}\) being the undeformed and current positions of a material point, and \(J=\det \varvec{{\varvec{F}}}\). The coefficients G and K are the shear and bulk modulus, respectively, connected via the Poisson ratio \(\nu \): \(K=\frac{2G}{3}\frac{1+\nu }{1-2\nu }\). Finally, the governing equation of solid deformation is given by:

$$\begin{aligned} \nabla \cdot \varvec{{\varvec{\sigma }}} = 0. \end{aligned}$$
(2)

In terms of boundary conditions, the basal surface of the basement membrane has zero displacement, while all the other surfaces are free to deform with no load or constraint (Fig. 2). In the paracellular geometry of Fig. 2a, we deploy distributed elastic springs between the opposite lateral surfaces of the two adjacent ECs to model the cell–cell adhesion force because of VE-cadherin bonds. Details of the junction model will be given below.

We use the augmented Lagrange method and the penalty method to model the normal contact force between the rod and the EC surface during transcellular and paracellular penetration, respectively [18]. Also we assume a frictionless contact between the protrusion and the ECs [19]. The sets of equations were solved using COMSOL Multiphysics. A special numerical challenge is to capture the very large strain under the ILP in the transcellular setup, and we meet it by employing a very high-aspect-ratio rectangular mesh in the surrounding area with quintic Lagrange shape function for the elements. Details can be found in the online Supplemental Information (SI).

2.3 Parameter estimation

To test the tenertaxis hypothesis, a prerequisite is to estimate the mechanical rigidity of the EC body as well as the EC junctions. These bear directly on the resistance against transmigration through the transcellular and paracellular routes, and hence the preference for one route or the other. In the following, we estimate the elastic modulus of the human lung microvascular endothelial cells and the mechanical properties of their junction. The outcome of this exercise is the baseline values summarized in Table 1 for some key mechanical parameters.

Table 1 Baseline values for key parameters used in our model, estimated for human lung microvascular endothelial cells
Full size table

2.3.1 Endothelial cells

The elastic modulus of ECs has most often been measured from the response of the cells to indentation, e.g., by the tip of an atomic force microscope (AFM) [24]. Generally, the elastic modulus tends to be higher above the nucleus and lower toward the periphery. For different types of human tissues, the EC body modulus, away from the nucleus, falls in the range \(E= 1000\sim 5800\) Pa [17, 20, 24,25,26,27]. For example, Viswanathan et al. [20] reported \(E=1800\) Pa near the periphery of human lung microvascular endothelial cells (HLMVECs). For human umbilical vein endothelial cells (HUVECs), E varies from 6800 Pa on the top of the nucleus to 1400 Pa near the cell edge [25]. For human aortic endothelial cells (HAECs), \(E=5800\) Pa is recorded above stress fibers while \(E=1500\) Pa away from stress fibers [27]. To propose the tenertaxis hypothesis, Martinelli et al. [3] quoted data on HLMVECs among other types of ECs. Accordingly, we adopt the HLMVEC cell body modulus [20] for simulating the paracellular protrusion: \(E_\mathrm{para} = 1800\) Pa. This value is also close to the modulus of the cell body of HUVECs and HAECs.

To estimate the EC modulus during transcellular diapedesis, a complication arises from the fact that a leukocyte, upon adhering to the apical surface of an EC, actively remodels the EC cortex underneath [21, 28]. Thus, the effective EC modulus is much reduced in comparison to areas away from the leukocyte attachment site. After removing an attached leukocyte by nano-surgery, Isac et al. [21] observed depolymerization of F-actin and remodeling of the EC cytoskeleton underneath the site of attachment. Barzilai et al. [28] showed similar disassembly of actin filaments underneath a leukocyte protrusion. Furthermore, Isac et al. [21] measured the EC elastic modulus at the attachment site using AFM indentation. The site can be up to 10 times softer than its surrounding area, with E dropping from \(\sim 3000\) to \(\sim 300\) Pa for HUVECs. This trend is consistent with the observation that indenting HAECs above stress fibers would yield a much stiffer modulus than away from stress fibers [27]. Unfortunately, we have found no quantitative data on the softening of HLMVEC modulus. Since we wish to model the HLMVEC-based experiments of Martinelli et al. [3], we borrow the softening factor of 10 from the HUVEC data [21]. Based on the modulus of 1800 Pa for the intact HLMVEC, we adopt an effective \(E_\mathrm{tran} = 180\) Pa for the transcellular simulation.

For both paracellular and transcellular simulations, the elastic modulus of the basement membrane is taken to be \(E_\mathrm{BM}=5000\) Pa, similar to the elastic modulus of extracellular matrix [17].

2.3.2 Endothelial junctions

The EC junction consists of a narrow gap bridged by molecular bonds (Fig. 1). As an ILP penetrates through the gap, it widens the gap by separating the EC surfaces and stretching the bonds. Following earlier studies, we represent the bonds by linear elastic springs [9, 10, 14], with the following elastic force on each bond:

$$\begin{aligned} f_{b}=k_{b} (L_{b}-L_{0} ), \end{aligned}$$
(3)

where \(k_{b}\) is the spring constant, \(L_{b}\) is the bond length and \(L_0\) is the equilibrium bond length. The spring constant \(k_{b}\) was reported to be \(k_b = 10^{-5}\sim 10^{-2}\) N/m for cell–cell adhesion bonds [9] and the equilibrium bond length \(L_0 = 10 \sim 20\) nm [6, 10]. We have taken \(k_b= 2\times 10^{-5}\) N/m and \(L_0=10\) nm to be our baseline values, the latter being consistent with the undeformed gap width at the junction.

The areal bond density \(N_b\) has been modeled by a kinetic equation for the formation and dissociation of bonds, with the rate constants being functions of tension in the bonds [10, 14]. For simplicity, we ignore the kinetics of \(N_b\) and adopt a constant equilibrium value based on prior observations. The total number of bonds between two ECs falls in a wide range, and the areal density of bonds can be estimated as \(N_b = 5 \sim 50{,}000\) \(\upmu \)m\(^{-2}\) [6, 9, 10]. Furthermore, from the measured normal stress (\(\sim \,1\) nN/\(\upmu \)m\(^2\)) at the junction [22] and the maximum force (\(\sim \) 50 pN) that each bond can sustain [23], one can estimate \(N_b = 20\) \(\upmu \)m\(^{-2}\). This is the value that we have adopted.

In our finite element simulation, the bond force is applied to each mesh point on the portion of the EC surface that forms the junction. Each mesh point is assigned an associated surface area \(A_i\) based on the mesh configuration. The bond length \(L_b\) at the mesh point is the horizontal distance to the opposite cell face. Thus, we calculate a horizontal bond force \(F_i= A_i N_b k_b (L_b-L_0)\) and apply it onto the mesh point.

For the parameters estimated above, the junctional springs turn out to contribute a small amount of the resistance to paracellular diapedesis; most of the resistance comes from the deformation of the EC cells. This is because the bond forces are horizontal, and only contribute indirectly to the resistance by affecting the EC surface shape. In particular, its contribution vanishes toward the end, when the junctional gap becomes widened more or less uniformly to the size of the rod. The minor role of the junctional springs alleviates the concern about uncertainties in estimating the above parameters.

We have carried out a parametric study to examine how the model predictions vary with the key geometric and mechanical parameters. Details are given in the online SI.

3 Results and discussion

We focus on the resistance force on the ILP as it extends across the endothelium through the transcellular or paracellular route. This resistance force is used to predict which is the preferred route with the smaller resistance. We first present results using the baseline parameters corresponding to wild-type HLMVECs and then investigate the effects of manipulating the elastic modulus of the endothelium.

3.1 Transmigration predicted using the baseline parameters

3.1.1 Paracellular route

In our quasi-static protocol, we move the rod representing the ILP downward into the junctional gap by 10 nm each time and compute the resistance force \(F_r\) on it from the contact surfaces between the rod and the two neighboring ECs. Figure 3a plots the force as a function of the displacement of the tip of the rod \(d_\mathrm{rod}\). Because of the geometry of the problem, the rod travels downward by about 75 nm before it makes contact with the EC surfaces at the top of the junctional gap. Initially, the rod deforms the upper edges of the ECs at the junction and opens it up into a wedge. The resistance force \(F_r\) rises sharply as a result. It peaks at 35 pN, for \(d_\mathrm{rod}\approx 0.25\) \(\upmu \)m, when the rod pushes the edges of the ECs apart so that the tip of the wedge reaches the basal EC surface (Fig. 3b). As the rod pushes further down, it opens an increasingly larger portion of the junction into a cylindrical hole, which does not provide much upward resistance. Therefore, \(F_r\) starts to decline with \(d_\mathrm{rod}\). This suggests a discontinuous jump in a dynamic simulation, and we will return to this point later. The minimum resistance \(F_r=20.3\) pN occurs at \(d_\mathrm{rod}\approx 0.54\) \(\upmu \)m (Fig. 3c). Afterward, the basal face of the EC starts to cause appreciable deformation in the basement membrane, which is much stiffer than the EC (Table 1). Thus, \(F_r\) starts to increase again with displacement. This increase continues till the end of the penetration, when the rod reaches the basement membrane (\(d_\mathrm{rod} = d_b=1\) \(\upmu \)m). Thus, the maximum resistance force \(F_\mathrm{max}= 45\) pN occurs at the end of the paracellular transmigration. The process is illustrated in Movie 1 in the SI.

Fig. 3
figure 3

a Vertical contact force as a function of the rod displacement during paracellular protrusion. The dashed arrow suggests how the rod will pass over an unstable portion of the curve in a dynamic simulation with a prescribed force. Penetration occurs when \(d_\mathrm{rod}\) reaches the breakthrough threshold \(d_b=1\) \(\upmu \)m. b Cross-sectional view of the mid-plane normal to the groove atop the EC junction (see Fig. 2a) at several points of the transmigration, with contours of the von Mises stress in Pa. The ECs and the basement membrane are colored dark blue and light blue, respectively

Full size image

Is it reasonable to expect the ILP to produce a protrusive force large enough to overcome such a resistance? Earlier studies showed that the protrusive force is mainly due to actin filaments polymerizing in the core of the ILP, with little direct contribution from myosin [29]. The force generated by a single polymerizing actin filament ranges between 1.3 and 9 pN [30, 31]. Each protrusion may contain between 10 and 30 actin filaments [32]. Thus, the ILP can easily generate a protrusive force that matches the resistance of the paracellular route.

Because of the quasi-static setup of the problem, in which we prescribe the displacement of the rod, the force–displacement curve of Fig. 3a has interesting implications for a “dynamic simulation” where we prescribe the pushing force \(F_r\) on the rod, and compute its displacement \(d_\mathrm{rod}(t)\). If \(F_r\) is below 35 pN, the rod will stop on the first upslope of the curve in Fig. 3a. With larger forces, the rod passes the peak and then jumps over an unstable portion of the curve to land on a larger displacement on the second upslope of the curve. This is indicated by the red dashed arrow in Fig. 3a. In fact, we have done such dynamic simulations to confirm the above scenarios. Details are discussed in the SI, with a simulation shown in Movie 2. Of course, there would be a jump in the opposite direction if we start with an initial state of full penetration and then gradually decrease the protrusion force. But such a scenario is not relevant to the diapedesis process.

3.1.2 Transcellular route

To test transcellular protrusion, we again adopt a quasi-static approach by moving the ILP downward by small increments and compute the elastic resistance by the EC. Figure 4a plots the vertical contact force during the progression of the rod, while Fig. 4b depicts the contours of the von Mises stress on the meridian plane at several points of the transcellular penetration of the EC. The protrusive force on the EC produces a transcellular tunnel similar to the fingerlike protrusions of leukocytes observed in the experiments [3, 4]. The force \(F_r\) increases monotonically but nonlinearly with the depth of penetration \(d_\mathrm{rod}\), owing to the geometric nonlinearity of the deformation and the hyperelastic behavior of the EC. Thus, when the protrusion reaches the threshold of \(d_b=0.98\) \(\upmu \)m, the resistance force attains its maximum during the transcellular penetration: \(F_\mathrm{max}=72\) pN.

Fig. 4
figure 4

a Vertical contact force as a function of the rod displacement during transcellular protrusion. The vertical dashed line marks the breakthrough threshold \(d_b=0.98\) \(\upmu \)m. b The meridian plane of the axisymmetric geometry (see Fig. 2b) at several points of the transmigration, with contours of the von Mises stress in Pa. The ECs and the basement membrane are colored dark blue and light blue, respectively

Full size image

Prior experiments show that paracellular transmigration is the preferred mode for HLMVECs [3, 4] and HUVECs [33]. Across an HLMVEC monolayer, for example, paracellular diapedesis occurs about 65% of the time in vitro [3, 4]. Based on baseline parameters corresponding to HLMVECs, our model predicts a maximum mechanical resistance of 45 pN and 72 pN in the paracellular and trancellular routes, respectively. This shows that the penetration is easier using the paracellular route, in agreement with experimental observations [3, 4]. To probe the tenertaxis hypothesis further, however, we need to manipulate the relative level of resistance in the EC body and junctional areas. This was achieved in prior experiments by drug or hormone treatment, shear flow effect and comparing different cell types that offer different levels of junctional resistance. In the following, we will simulate such changes by manipulating the endothelial resistances at the cell junctions.

3.2 Effect of manipulating endothelial stiffness at the junction

In a series of experiments, Martinelli et al. [3, 4, 34] demonstrated two interesting trends. First, for HLMVECs in their natural physiologic state, the paracellular route is preferred and accounts for 65% of the total transmigration events. Second, by strengthening the EC stiffness near junctions, a reversal in this preference can be achieved, with a majority (70%) of the transmigration occurring through the transcellular route. This has motivated us to manipulate the junctional strength in our model to see how that modifies the relative resistance of the two routes.

Using the hormone adrenomedullin, Martinelli et al. [3] managed to increase the level of cortical F-actin near the HLMVEC edges. As a result, the EC elastic modulus rose to 1.75 times of its baseline value. Moreover, Viswanathan et al. [20] reported increase of HLMVEC modulus by a factor of 2.33 after treatment with hepatocyte growth factor. Away from the cell edge, the EC body modulus is not affected by the chemical treatment, and neither is the transcellular migration [3, 20, 35, 36]. In view of these experimental data, we will increase \(E_\mathrm{para}\) by a factor of 1.75 and 2.33 from its baseline value of 1800 Pa, to \(E_\mathrm{para}=3150\) Pa and 4200 Pa, respectively. None of the experiments indicated modifications of the VE-cadherin-based molecular bonds at the EC junctions. Thus, we have kept the baseline values of \(k_{b}=2\times 10^{-5}\) N/m.

Fig. 5
figure 5

Effect of raising the EC modulus near the cell edges on the resistance force during paracellular penetration. The baseline value is \(E_\mathrm{para}=1800\) Pa (see Table 1), and the two higher values correspond to experimental manipulations of EC edge rigidity [3, 20]

Full size image

Figure 5 depicts the effect of stiffer EC-edge modulus on the resistance force during ILP penetration. The maximum resistance force \(F_\mathrm{max}\) has increased markedly for the strengthened junctions, from \(F_\mathrm{max}=45\) pN for the wild-type EC to 74 pN for \(E_\mathrm{para}=3150\) Pa and further to 90 pN for \(E_\mathrm{para}=4200\) Pa. For both cases of elevated \(E_\mathrm{para}\), the resistance of paracellular transmigration now surpasses that of the transcellular route. Thus, the model reproduces the observations of Martinelli et al. [3] that stiffening the EC near cell junctions can make the transcellular route preferable. Our mechanical tests have tentatively confirmed the hypothesis of tenertaxis.

4 Conclusion

As part of the immune reaction, leukocytes transmigrate across the endothelium either directly through the body of an endothelial cell (EC) (transcellular route) or through the junction between adjacent ECs (paracellular route). To rationalize the preference between the two routes, Martinelli et al. [3] hypothesized that leukocytes seek the path of least mechanical resistance in a process called tenertaxis. In this study, we have examined this hypothesis using numerical computation of the mechanical resistance encountered by the leukocyte protrusion during paracellular or transcellular penetration. Model predictions show that normally the paracellular route is the preferred route for human lung endothelia. Using model parameters corresponding to the human lung microvasculature, our computations show that the required force to penetrate the endothelium in the transcellular route (\(F_\mathrm{max}=72\) pN) is greater than that of the paracellular route (\(F_\mathrm{max}=45\) pN). This rationalizes the preference of leukocytes to use the paracellular route most of the time.

Motivated by experiments that enhanced the junctional integrity of endothelium through the use of modifying agents, we have demonstrated that by increasing the elastic modulus of the EC near the junction, the mechanical resistance of the paracellular route may surpass that of the transcellular route. This will make the transcellular route preferable, in agreement with experimental observations [3]. Thus, our mechanical tests have tentatively confirmed the hypothesis of tenertaxis.

Our model is purely mechanical and aims at testing the mechanical feasibility of the tenertaxis hypothesis. Thus, it has incorporated many assumptions and simplifications. First, we have ignored all biochemical signaling in the complex process of diapedesis. In particular, we have neglected the kinetics of F-actin polymerization, which determines the force generation inside the invadosome-like protrusion (ILP), and the breakage of molecular bonds in the EC junction during paracellular transmigration. A more complete model should integrate such biochemical kinetics with the mechanical deformation of the EC. Second, our model does not account for the leukocyte nucleus. Instead, we follow the experimental work [2,3,4] in focusing solely on the ILP as the “probe” for measuring the resistance of various endothelial components. Because of the complex morphology and dynamic behavior of the leukocyte nucleus, its role in diapedesis is an open and actively debated question [37, 38]. Third, we model the EC as hyperelastic and disregard the cytoplasm and potential viscoelastic responses of the endothelium [39]. Finally, we have treated the ILP as an inert solid object. Experimental evidence suggests that there is a dynamic interaction between the leukocyte and the ECs. Through mechanical and chemical signaling pathways, the leukocyte may induce cortical contraction inside the EC to open up the cell junctions and facilitate paracellular passage [40, 41]. In return, the EC appears to modulate the lifetime of the ILP through an active feedback mechanism, such that unsuccessful protrusions retract quickly while more successful ones grow and persist much longer [33]. Thus, it is clear that leukocyte diapedesis is a complex process that involves a rich array of biochemical and mechanical processes. Tenertaxis is a simple but promising idea that explains the outcome from predominantly mechanical factors. Our modeling offers quantitative support for this concept.

We end by noting the clinical relevance of the above discussion. Endothelial barrier function is perturbed in many disease states (e.g., sepsis, COVID-19 and atherosclerosis, to name a few) that involve endothelial activation and enhanced leukocyte transit through the paracellular space. Once in the tissue, the leukocytes activate and promote inflammation and cause organ damage. The use of glucocorticoids, which “stiffen” the endothelial cell junctions by enhancing synthesis of junction proteins such as occludin and VE-cadherin, reduces leukocyte entry and dampens inflammation [42]. Endothelial barrier function is thus an important target of therapeutic discoveries for many inflammatory conditions.