1. Introduction
Interfacial flows with insoluble or soluble surfactants are intensely studied in the recent literature and the presence of the surface-active agent has a decisive role in the stability characteristics and transitions, as well as in accompanying transport phenomena (Frenkel & Halpern Reference Frenkel and Halpern2002; Blyth & Pozrikidis Reference Blyth and Pozrikidis2004; Pereira et al. Reference Pereira, Trevelyan, Thiele and Kalliadasis2007; Craster & Matar Reference Craster and Matar2009; Kalogirou & Papageorgiou Reference Kalogirou and Papageorgiou2015; Manikantan & Squires Reference Manikantan and Squires2020; Hu, Fu & Yang Reference Hu, Fu and Yang2020; Constante-Amores et al. Reference Constante-Amores, Batchvarov, Kahouadji, Shin, Chergui, Juric and Matar2021; D’Alessio & Pascal Reference D’Alessio and Pascal2021; Samanta Reference Samanta2021; Miguet, Dorbolo & Scheid Reference Miguet, Dorbolo and Scheid2024; Samanta Reference Samanta2025). The dynamics of such flows is frequently dictated by the coupling of the flow field with the characteristics of the surfactant monolayer, for example its surface elasticity and viscosity, and with the kinetics of surfactant exchange between the surface and the bulk. This coupling becomes more involved when the monolayer is dense, or even close to collapse, in which case its dynamics cannot be described by simple linear relations.
A relevant application, of particular interest to mammalian physiology, is the behaviour of the surfactant-laden liquid layer that lines internally the lung airways and alveoli (Gaver & Grotberg Reference Gaver and Grotberg1990; Halpern, Jensen & Grotberg Reference Halpern, Jensen and Grotberg1998; Matar, Zhang & Craster Reference Matar, Zhang and Craster2003; Grotberg Reference Grotberg2011; Kim et al. Reference Kim, Choi, Zasadzinski and Squires2011; Filoche, Tai & Grotberg Reference Filoche, Tai and Grotberg2015; Muradoglu et al. Reference Muradoglu, Romanò, Fujioka and Grotberg2019; Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022; Romanò et al. Reference Romanò, Muradoglu and Grotberg2022). In the case of alveoli, this layer consists of an aqueous phase, with estimated thickness
$0.1{-}1 \; {\unicode{x03BC}\textrm {m}}$
depending on lung inflation and health condition (Bastacky et al. Reference Bastacky, Lee, Goerke, Koushafar, Yager, Kenaga, Speed, Chen and Clements1995; Wei et al. Reference Wei, Fujioka, Hirschl and Grotberg2005). Its interface with the air is covered by a monolayer of pulmonary surfactant, a mixture of lipids (mainly dipalmitoylphosphatidylcholine (DPPC)) and proteins (SP-A, SP-B, SP-C Zuo et al. (Reference Zuo, Veldhuizen, Neumann, Petersen and Possmayer2008)). These molecules are sparingly soluble, so any excess resides inside the liquid film in the form of micelles (Perez-Gil Reference Perez-Gil2008). The key property of the surfactant monolayer is its ability to reduce drastically surface tension at the interface, thus making the alveoli more compliant and reducing the metabolic work of breathing (Zasadzinski et al. Reference Zasadzinski, Ding, Warriner, Bringezu and Waring2001; Rugonyi, Biswas & Hall Reference Rugonyi, Biswas and Hall2008; Zhang et al. Reference Zhang, Wang, Fan and Zuo2011).
The hydrodynamics of the liquid layer lining the alveoli has been repeatedly a topic of investigation during recent decades (Gradon & Podgorski Reference Gradon and Podgorski1989; Podgorski & Gradon Reference Podgorski and Gradon1993; Espinosa & Kamm Reference Espinosa and Kamm1997; Zelig & Haber Reference Zelig and Haber2002; Wei et al. Reference Wei, Benintendi, Halpern and Grotberg2003, Reference Wei, Fujioka, Hirschl and Grotberg2005; Halpern et al. Reference Halpern, Fujioka, Takayama and Grotberg2008; Kang et al. Reference Kang, Chugunova, Nadim, Waring and Walther2018; Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022). The reason for this interest is that slow shearing motions (i.e. velocities parallel to the alveolar epithelium), which may develop triggered by the radial oscillation of the alveolar wall, are potentially of importance for lung homeostasis, for example in cell–cell signalling and in cleansing the alveolus from deposited particles. Two key mechanisms that may create a shearing motion are Marangoni (elastic) stresses that result from spatial variation of surface tension and capillary stresses that result from spatial variation of the interfacial curvature.
Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) modelled alveoli as spherical caps in self-similar oscillations and argued that the boundary conditions imposed at the rim are critical for correct modelling of the flow field in the liquid lining. A novel boundary condition was implemented by these authors, which takes into account the small but finite thickness of the mid-alveolar septum and avoids the imposition of pre-specified values of film thickness and surfactant concentration at the rim. This approach appears to provide realistic estimates for key variables, such as the shear stress imposed along the epithelium, and also offers interesting predictions for conditions in a diseased lung (Livanos et al. Reference Livanos, Bouchoris, Aslani, Gourgoulianis and Bontozoglou2024).
However, the study by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) is limited to linear and weakly nonlinear oscillations of the alveolar wall around its equilibrium position, i.e. formally to alveolar radius variations of order
$\,\pm 1\,\%$
. Thus, the following questions naturally arise: (i) Is linear theory sufficient for the description of smooth breathing, where radius oscillations are typically around
$3\,\%{-}5\,\%$
? (ii) Are there any significant nonlinear phenomena in this breathing regime? (iii) What is the flow field in the liquid layer during deep breathing, where the alveolar radius may vary by
$10\,\%{-}20\,\%$
? The investigation of such questions necessitates solution of the full nonlinear problem, which is the goal of the present work.
The pulmonary surfactant monolayer is known to undergo significant structural changes with the occurrence of large variations in its surface area. These changes have been extensively investigated in the relevant literature and many illuminating reviews are available (Possmayer et al. Reference Possmayer, Nag, Rodriguez, Qanbar and Schürch2001; Zasadzinski et al. Reference Zasadzinski, Ding, Warriner, Bringezu and Waring2001; Wüstneck et al. Reference Wüstneck, Perez-Gil, Wüstneck, Cruz, Fainerman and Pison2005; Lee Reference Lee2008). In particular, the monolayer may sustain large compressions during contraction, eventually reducing its surface tension to values below
$0.005\, \textrm {N m}^{-1}$
, while changing its structure to a liquid expanded phase and then to a tilted condensed phase (Crane, Putz & Hall Reference Crane, Putz and Hall1999; Yan, Piknova & Hall Reference Yan, Piknova and Hall2005; Zuo et al. Reference Zuo, Veldhuizen, Neumann, Petersen and Possmayer2008; Zhang et al. Reference Zhang, Wang, Fan and Zuo2011; Parra & Perez-Gil Reference Parra and Perez-Gil2015).
Further compression results in collapse of the surfactant monolayer and formation of a reservoir sublayer (Krueger & Gaver Reference Krueger and Gaver2000). The onset of collapse is characterised by a minimum value of the surface tension,
$\sigma _{min}$
, and a maximum value of the surface concentration of surfactant,
$\varGamma _{max}$
, both of which remain constant until the end of the contraction stage. The reservoir sublayer appears to form by buckling and folding events, giving rise to protrusions that extend away from the monolayer and towards the interior of the liquid while remaining connected to the monolayer (Lipp et al. Reference Lipp, Lee, Takamoto, Zasadzinski and Waring1998; Schief et al. Reference Schief, Antia, Discher, Hall and Vogel2003; Baoukina et al. Reference Baoukina, Monticelli, Risselada, Marrink and Tieleman2008; Lee Reference Lee2008; Saad, Neumann & Acosta Reference Saad, Neumann and Acosta2010; Phan, Lee & Shin Reference Phan, Lee and Shin2016; Xu, Yang & Zuo Reference Xu, Yang and Zuo2020). The connection of the reservoir to the monolayer is very important because it allows rapid replenishment of the monolayer content during expansion, thus restricting the rise in surface tension (Wüstneck et al. Reference Wüstneck, Wüstneck, Fainerman, Miller and Pison2001; Parra & Perez-Gil Reference Parra and Perez-Gil2015).
In the present work, the evolution equations for the liquid film thickness,
$h$
, and the surface surfactant concentration,
$\varGamma$
, of an internally lined oscillating spherical cap – derived by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) using a lubrication approximation in the Stokes limit – are supplemented by an evolution equation for the pure monolayer concentration,
$\varGamma _m$
. Thus, capillary and Marangoni stresses, which actually vary with the monolayer and not the total surface concentration, are correctly described for large-amplitude oscillations that result in monolayer collapse.
A novel kinetic expression is introduced, which models the replenishment of the monolayer by surfactant from the reservoir,
$ (\varGamma -\varGamma _m)$
, during re-expansion of the alveolus. Even more important, the kinetic expression includes a term that restricts monolayer concentration to values
$\varGamma _m\leqslant \varGamma _{max}$
, while introducing a parameter,
$k_{rc}$
, that quantifies the abruptness of monolayer collapse. The resulting system of nonlinear equations is spatially discretised by a Galerkin finite-element method and solved by a Newton–Raphson iterative scheme, with the time derivatives treated fully implicitly.
The paper is organised as follows. The problem set-up in § 2 includes formulation of the dynamics of the liquid layer, the dynamics and the mass balances of the surfactant and the boundary and initial conditions. In § 3, the equations are scaled using equilibrium variables and their numerical solution is outlined. Section 4 contains results for the dynamics of large-amplitude oscillations (without and with monolayer collapse) and investigates in particular the displacement of insoluble particles (passive tracers) deposited on the liquid layer and the drift of surfactant from the alveolar opening. Finally, conclusions are drawn in § 5.
2. Problem set-up
2.1. Outline of the modelling procedure
Model development necessitates combined consideration of (i) the fluid dynamics of the liquid layer, (ii) the mass exchange dynamics at the interface and (iii) the conservation equations of surfactant. The derivations, although not particularly difficult, are involved and with rather numerous intermediate steps, so the following outline may be helpful.
The fluid dynamics of the liquid layer is described by the Navier–Stokes equation in the low Reynolds number limit, with free-surface boundary conditions applied along the interface. The dynamic boundary condition contains the local values of surface tension (in the capillary stress term) and its spatial gradients (in the Marangoni stress term). The surface tension is considered a sole function of the local monolayer concentration,
$\varGamma _m$
, and thus its values – as well as those of surface elasticity – are derived from an equation of state that follows from Gibbs’ theory, with key novelty the inclusion of the intrinsic compressibility of the adsorbed monolayer.
Surfactant concentration on the interface is described by two terms,
$\varGamma _m$
, which gives the surface concentration of the adsorbed monolayer and
$\varGamma$
, which, includes both the monolayer surfactant and the surfactant contained in the sublayer reservoir that forms during monolayer collapse. Two conservation equations are thus formulated to follow the spatio-temporal evolution of the above surface concentrations, which include the contributions of interfacial expansion/contraction, of convection and diffusion along the interface and of mass exchange with the environment.
The mass exchange terms are crucial inputs to the surfactant conservation equations, and they consist of the following two contributions. The flux
$j_b$
, which describes adsorption/desorption between the monolayer and the bulk, and is modelled by Langmuir kinetics in terms of an adsorption rate and an equilibrium constant. The flux,
$j_r$
, which describes internal exchange between the monolayer and the reservoir, and is modelled by a Langmuir-like term that accounts for monolayer replenishment from the reservoir during interface expansion and by a singular term that limits maximum monolayer concentration below
$\varGamma _{max}$
during interface contraction. In the following paragraphs, each of the aforementioned dynamic components of the model is considered in detail.
2.2. Dynamics of the liquid layer
The alveolus is modelled as a spherical cap with an opening angle
$2 \theta _0$
(
$\theta _0=30^\circ$
for a normal alveolus), as shown in figure 1(a). During oscillation, the opening angle remains constant, i.e. the alveolus retains a self-similar shape (Tsuda, Henry & Butler Reference Tsuda, Henry and Butler2008). Although there are alternative choices for the exact alveolar shape, the truncated sphere remains the most common model geometry in the literature (Kolanjiyil & Kleinstreuer Reference Kolanjiyil and Kleinstreuer2019). The alveolar rim is taken to have a semicircular cross-section of small but finite radius
$r_0$
(figure 1
b). Thus, the local septal thickness is
${\sim}2r_0$
. Due to the finite thickness, water and surfactant accumulate on the rim and local values of film thickness and surfactant concentration, which continuously vary during alveolar oscillations, enter the boundary conditions (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022).
(a) A sketch of the spherical cap with the main problem parameters and (b) magnification of the rim with the local parameter values (red colour marks the hydrophobic end of surfactant molecules).
The liquid layer is Newtonian (Grotberg Reference Grotberg2011), with density
$\rho$
and viscosity
$\mu$
. The flow is taken to be symmetric in the circumferential direction and is analysed in a spherical coordinate system
$(r,\theta,\phi)$
located at the centre of the cap. Thus, the flow field is described by the velocity vector
$\boldsymbol{u}=(u_{r}(r,\theta,t),u_{\theta }(r,\theta,t),0)$
and the pressure distribution
$p(r,\theta,t)$
. The air–liquid interface is located at
$r=R(t)-h(\theta,t)$
, where
$h(\theta,t)$
is the liquid film thickness. It is further characterised by its surface tension,
$\sigma (\theta,t)$
, which varies in space and time as it is determined by the local concentration of the surfactant monolayer.
Following standard practice in the literature (Podgorski & Gradon Reference Podgorski and Gradon1993; Haber et al. Reference Haber, Butler, Brenner, Emmanuel and Tsuda2000; Zelig & Haber Reference Zelig and Haber2002; Wei et al. Reference Wei, Fujioka, Hirschl and Grotberg2005; Kang et al. Reference Kang, Chugunova, Nadim, Waring and Walther2018), the flow is described by the continuity and the quasi-steady Stokes equation, subject to the no-slip boundary condition on the wall and the kinematic and dynamic boundary conditions along the interface. The dynamic condition includes components normal and tangential to the interface, the former proportional to surface tension and local curvature (capillary stress) and the latter proportional to the local spatial gradient of surface tension (Marangoni stress). The equations are further simplified by a lubrication assumption, which is based on the large scale difference between the
$r\hbox{-}$
and
$\theta \hbox{-}$
directions (
$h\lt \lt R$
). The final result is the following evolution equation that describes the spatio-temporal variation of the liquid film thickness, driven by capillary and Marangoni stresses (respectively, the first and second terms in the parenthesis). The details of the mathematical derivation are presented in a previous paper (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022) and are also discussed by Kang et al. (Reference Kang, Chugunova, Nadim, Waring and Walther2018)
\begin{equation} \frac {\partial h}{\partial t} + \frac {2h\dot {R}}{R} - \frac {1}{R^2\sin \theta }\frac {\partial }{\partial \theta }\left (\frac {h^3\sin \theta }{3\mu }\frac {\partial p}{\partial \theta }-\frac {h^2\sin \theta }{2\mu }\frac {\partial \sigma }{\partial \theta }\right) = 0. \end{equation}
It is noted that, in the lubrication approximation, pressure is uniform in the radial direction and can be described by the following expression, using the normal boundary condition on the interface:
\begin{equation} p(\theta,t) = p_{air} - \sigma \left [\frac {2}{R}+\frac {2\,h}{R^2}+\frac {1}{R^2\sin \theta }\frac {\partial }{\partial \theta }\left (\sin \theta \,\frac {\partial h}{\partial \theta }\right)\right]. \end{equation}
The term in parenthesis in (2.1) is equal to
$-Q(\theta,t)/2\pi$
, where
\begin{equation} Q(\theta,t)=2\pi R\,\sin \theta \int _{R-h}^{R} u_{\theta }\,{\rm d}r=2\pi \sin \theta \left (-\frac {h^3}{3\mu }\frac {\partial p}{\partial \theta }+\frac {h^2}{2\mu }\frac {\partial \sigma }{\partial \theta }\right)\,, \end{equation}
is, in the lubrication limit, the volumetric flow rate along the entire
$\phi \hbox{-}$
circumference. Thus, (2.1) may be rewritten as
\begin{equation} \frac {\partial }{\partial t}\big (2\pi R^2\sin \theta \,h\big) + \frac {\partial Q}{\partial \theta } = 0, \end{equation}
which is intuitively more evident as the mass balance of liquid in a differential control volume along the entire
$\phi \hbox{-}$
periphery at the location
$\theta$
.
The derivations leading to (2.1) – as well as the entire paper – neglect the effect of the dynamics of airflow inside the alveolus during the breathing cycle. The following simple scaling argument, based on the continuity of shear stress and velocity at the air–liquid interface, shows that the effect is indeed negligible. More specifically, if
$u_{s,a}$
is the liquid interfacial velocity tentatively imposed by gas shear
\begin{equation} \mu _a\left (\frac {\partial u}{\partial y}\right)_a\sim \mu \frac {\partial u}{\partial y}\Rightarrow \mu _\alpha \frac {u_a-u_{s,a}}{\bar {R}}\sim \mu \frac {u_{s,a}-0}{\bar {H}}\Rightarrow u_{s,a}\sim u_a\frac {\mu _a}{\mu }\frac {\bar {H}}{\bar {R}}. \end{equation}
Given that
$\mu _a/\mu \sim 1:50$
and
$\bar {H}/\bar {R}\sim 0.3:100$
, the liquid velocity imposed by air shear is four orders of magnitude smaller than the typical air velocity inside the alveolus, i.e. two to three orders of magnitude smaller than the shearing velocity on the liquid film triggered by Marangoni stresses (as shown by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022), the shearing velocity is one to two orders of magnitude smaller than the radial velocity of the alveolar wall).
An assumption, implicit in the above formulation, is the neglect of the potential role of surfactant aggregates in the liquid film dynamics. The dynamics of the aggregates themselves is probably related to the kinematics of surfactant exchange with the interface and is not in direct interest to the present study. However, given the relative sizes of the aggregates (
${\sim} 100\,\unicode{x03BC}\rm m$
) and the liquid film (
${\sim} 200{-}1000\,\unicode{x03BC}\rm m$
), it is reasonable to expect that, with increasing number density, the aggregates will have a non-negligible effect on the liquid dynamics. Possibilities for including this effect in more refined simulations are mentioned in the Concluding Remarks.
2.3. Dynamics of the surfactant
The pulmonary surfactant is a mixture of lipids and proteins, which are practically insoluble in water. It is indicative that the critical micelle concentration (CMC) of the major phospholipids (i.e. the maximum concentration of monomer beyond which monomers combine to form micelles) is
$10^{-9}{-}10^{-10}\,\textrm{M}$
(Marsh & King Reference Marsh and King1986). As a result, a portion of the surfactant populates the interface with local surface concentration
$\varGamma (\theta,t)$
, while the excess amount resides in the bulk in the form of aggregates (Zuo et al. Reference Zuo, Veldhuizen, Neumann, Petersen and Possmayer2008; Bykov et al. Reference Bykov, Milyaeva, Isakov, Michailov, Loglio, Miller and Noskov2021). There is ample experimental evidence (Baoukina & Tieleman Reference Baoukina and Tieleman2011; Parra & Perez-Gil Reference Parra and Perez-Gil2015; Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2021) indicating that the adsorption/desorption of surfactant to/from the interface occurs through its direct interaction with these aggregates. The pulmonary surfactant mixture is currently modelled by a single generic surfactant, which mimics the dynamic behaviour of the actual preparation (Wüstneck et al. Reference Wüstneck, Perez-Gil, Wüstneck, Cruz, Fainerman and Pison2005).
The description of the interfacial surfactant composition becomes more involved when the alveolar oscillations are strong enough to include the stages of monolayer collapse and subsequent replenishment. The surface concentration,
$\varGamma$
, now consists of the sum of the monolayer and reservoir concentrations,
$\varGamma =\varGamma _{m}+\varGamma _{r}$
. During compression, the monolayer concentration does not increase beyond the maximum value that corresponds to the minimum attainable surface tension. Thus,
$\varGamma _{m}\rightarrow \varGamma _{max}$
, where
$\sigma (\varGamma _{max})=\sigma _{min}$
, and any further increase of
$\varGamma$
results in an equal increase of
$\varGamma _{r}$
. It is important to re-iterate that, capillary and Marangoni stresses vary with
$\varGamma _m$
and not with
$\varGamma$
. Figure 2 provides a sketch of the interfacial region and the surfactant exchange processes during compression (figure 2
a) and during expansion (figure 2
b).
Schematic of the interface with the monolayer and the reservoir protrusions, and the bulk with surfactant aggregates. (a) Compression that leads to monolayer collapse and growth of the reservoir. (b) Re-expansion, when the monolayer concentration falls below
$\varGamma _{max}$
and is replenished by surfactant from the reservoir. Red colour marks the hydrophobic end of surfactant molecules and bars in (b) (not to scale) indicate the sizes of the monomer, the surface layer and the aggregates as measured by Xu et al. (Reference Xu, Yang and Zuo2020).
The local value of surface tension,
$\sigma$
, along the interface is a function of the local monolayer concentration, i.e.
$\sigma =\sigma (\varGamma _m)$
, and is described by the equation of state to be developed shortly. Considering a space coordinate,
$s$
, along the interface, the variation of surface tension and monolayer concentration are related by the chain rule, i.e.
$\partial \sigma /\partial s=({\rm d}\sigma /{\rm d}\varGamma _m)(\partial \varGamma _m/\partial s)$
. The sensitivity of
$\sigma$
to
$\varGamma _m$
is expressed by the Gibbs elasticity,
$E$
, where
\begin{equation} E = -\frac {{\rm d}\sigma }{{\rm d}\ln \varGamma _m} = -\varGamma _m \frac {{\rm d}\sigma }{{\rm d}\varGamma _m}. \end{equation}
A novel feature of the surfactant model is the inclusion of an intrinsic compressibility,
$\alpha$
, of the adsorbed monolayer, as defined and justified by Fainerman, Miller & Kovalchuk (Reference Fainerman, Miller and Kovalchuk2002) and Kovalchuk et al. (Reference Kovalchuk, Loglio, Fainerman and Miller2004, Reference Kovalchuk, Miller, Fainerman and Loglio2005). More specifically, the molar surface area of the monolayer,
$\varOmega$
, is taken to vary linearly with surface pressure,
$\varPi =\sigma _0-\sigma$
, according to the relation
\begin{equation} \varOmega = \varOmega _{0} (1-\alpha \varPi), \end{equation}
where
$\sigma _0$
is the surface tension of pure water and
$\varOmega _0$
is the molar area at zero surface pressure. This correction is particularly significant for dense monolayers and was recently shown to offer quantitative agreement of the model dynamics with laboratory measurements using actual pulmonary preparations (Saad et al. Reference Saad, Neumann and Acosta2010; Xu et al. Reference Xu, Yang and Zuo2020).
In terms of
$\varOmega$
, the monolayer coverage,
$\gamma$
, is
$\gamma =\varGamma _m \varOmega$
, and thus the surface concentration at interfacial saturation,
$\varGamma _{m,\infty }$
(which corresponds to
$\gamma =1$
), varies with surface pressure, and is given by
\begin{equation} \varGamma _{m,\infty } = \frac {1}{\varOmega } = \frac {1}{\varOmega _{0}(1-\alpha \varPi)}. \end{equation}
Using Gibbs’ theory, which is valid for a Gibbs dividing surface and an ideal bulk phase, the following equation of state is derived (Kovalchuk et al. Reference Kovalchuk, Loglio, Fainerman and Miller2004), connecting surface pressure to interfacial coverage:
\begin{equation} \frac {\varPi \varOmega _0}{\mathcal{R}\mathcal{T}}\left (1-\alpha \frac {\varPi }{2}\right) = -\ln (1-\gamma),\end{equation}
where
$\mathcal{R}$
is the universal gas constant and
$\mathcal{T}$
the absolute temperature.
The surface concentration of surfactant varies continuously during breathing as a result of the periodic expansion and contraction of the interface. The variation of
$\varGamma _m$
from the equilibrium value,
$\varGamma _{eq}$
(corresponding to a non-oscillating alveolus for which
$\varGamma _{eq}\equiv \varGamma _{m,eq}$
) causes a mass flux,
$j_b$
, of surfactant between the monolayer and the bulk. Assuming a Langmuir adsorption isotherm, equilibrium is given by the expression
\begin{equation} KC_{10}=\frac {\varGamma _{eq}}{\varGamma _{eq,\infty }-\varGamma _{eq}},\end{equation}
where
$\varGamma _{eq,\infty }$
is the saturation value at equilibrium conditions,
$K$
is the equilibrium constant and
$C_{10}$
the CMC of the monomer in the bulk. Thus, mass exchange with the bulk is expressed as follows in terms of an adsorption rate,
$k_{ads}$
:
\begin{equation} j_b = k_{ads}C_{10}\left [(\varGamma _{m,\infty }-\varGamma _m)-\frac {\varGamma _m}{K\,C_{10}} \right], \end{equation}
while for an insoluble surfactant
$j_b=0$
. It is noted that the above expression for mass exchange is based on kinetic resistance at the interface (rather than diffusion in the bulk), an assumption strongly supported by the literature (Ingenito et al. Reference Ingenito, Mark, Morris, Espinosa, Kamm and Johnson1999; Saad et al. Reference Saad, Neumann and Acosta2010).
Mass exchange between the monolayer and the reservoir,
$j_r$
, changes direction during the breathing cycle, because the reservoir stores surfactant during compression (so that
$\varGamma _m$
does not exceed
$\varGamma _{max}$
) and releases it back to the monolayer during re-expansion of the alveolus (when
$\varGamma _m$
falls below
$\varGamma _{max}$
). This dual role is described by the following expression:
\begin{equation} j_{r} = k_r ( \varGamma _{m,\infty } - \varGamma _{m})\left (\varGamma - \varGamma _m\right) - \frac {k_{rc}}{\left (\varGamma _{max}-\varGamma _m\right)^2}. \end{equation}
The first term in (2.12) gives the rate of replenishment of the monolayer during expansion and is proportional to the potential for replenishment,
$ ( \varGamma _{m,\infty } - \varGamma _{m})$
, and to the available surfactant concentration in the local reservoir,
$ (\varGamma -\varGamma _m)$
. The second term limits the monolayer concentration to be below the maximum value
$\varGamma _{max}$
during compression. Parameter
$k_{rc}$
may be viewed as a measure of the stiffness of the dense monolayer, in the sense that the smaller
$k_{rc}$
the closer the monolayer concentration is to its limiting value without collapsing and the more abrupt the collapse. By choosing a very small value
$k_{rc}=10^{-5}$
, the contribution of the second term in (2.12) becomes significant only when
$\varGamma _m$
approaches close enough to
$\varGamma _{max}$
.
2.4. Mass balances of surfactant
The aforementioned variables,
$\varGamma$
and
$\varGamma _m$
, are functions of location and time, and their evolution is determined by appropriate mass conservation equations. Such equations must include convection and diffusion along the interface and mass exchange with the bulk and the reservoir, and must also take into account the effect of expansion/contraction of the interfacial area on surface concentration. It is noted that there is no need for a mass balance in the bulk, because the typical surfactant loading is much higher than the CMC of the monomer. Thus, the bulk is always saturated with monomer and the effect of bulk diffusion is negligible.
First, the case of weak alveolar oscillations – corresponding to shallow breathing manoeuvres – is considered. In this case, the surface concentration never approaches
$\varGamma _{max}$
, so
$\varGamma _m\equiv \varGamma$
. Thus, a single conservation equation is required, which – when simplified by a lubrication approximation – leads to the following result (Kang et al. Reference Kang, Chugunova, Nadim, Waring and Walther2018; Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022):
\begin{equation} \begin{aligned} &\frac {\partial \varGamma }{\partial t} + \frac {2\varGamma \dot {R}}{R} - \frac {1}{R^2\sin \theta }\frac {\partial }{\partial \theta }\left (\frac {\varGamma h^2\sin \theta }{2\mu }\frac {\partial p}{\partial \theta }-\frac {\varGamma h\sin \theta }{\mu }\frac {\partial \sigma }{\partial \theta }\right)\\ &\qquad - \frac {\mathcal{D}_s}{R^2\sin \theta }\frac {\partial }{\partial \theta }\left (\sin \theta \,\frac {\partial \varGamma }{\partial \theta }\right) = j_b. \end{aligned} \end{equation}
Here,
$\mathcal{D}_s$
is the surface diffusivity of surfactant and the terms in the two parentheses represent, respectively, the contributions of convection and surface diffusion. Eqaution (2.13) may take the following condensed form:
\begin{equation} \frac {\partial }{\partial t} \big(2\pi R^2\sin \theta \;\varGamma \big) + \frac {\partial Q_{\varGamma }}{\partial \theta } = 2\pi R^2 \sin \theta \, j_b, \end{equation}
where
$Q_{\varGamma }(\theta,t)$
is the interfacial flow rate of surfactant along the entire
$\phi -$
circumference
\begin{equation} Q_{\varGamma }(\theta,t)=2\pi R\,\sin \theta \left (u_s\varGamma -\frac {\mathcal{D}_s}{R} \frac {\partial \varGamma }{\partial \theta }\right),\end{equation}
and
$u_s$
is the interfacial water velocity
\begin{equation} u_s\equiv u_{\theta }(R-h,\theta,t)=-\frac {1}{2\mu }\frac {\partial p}{\partial \theta }\frac {h^2}{R}+\frac {1}{\mu }\frac {\partial \sigma }{\partial \theta }\frac {h}{R}. \end{equation}
In order to describe deep-breathing manoeuvres, i.e. alveolar oscillations that include monolayer collapse and replenishment, the following modifications are necessary. As noted previously,
$\sigma =\sigma (\varGamma _m)$
and
$E=E(\varGamma _m)$
, i.e. capillary and Marangoni forcing are both driven by the monolayer concentration and not the total surface value,
$\varGamma$
. Thus, convective and diffusive fluxes have as a driving force the monolayer concentration
$\varGamma _m$
. However, the reservoir extensions are attached to the monolayer (Parra & Perez-Gil Reference Parra and Perez-Gil2015), and consequently are transported with it. Thus, for each mol of monolayer surfactant transported by convection or diffusion
$\varGamma _r/\varGamma _m$
also follow, or equivalently a total of
$(\varGamma _r+\varGamma _m)/\varGamma _m=\varGamma /\varGamma _m$
is transported. As a result, the interfacial flow rate,
$Q_{\varGamma }$
, now becomes
\begin{equation} \begin{aligned} Q_{\varGamma }(\theta,t)=2\pi R\,\sin \theta \left [u_s(\varGamma _m)\,\varGamma _m-\frac {\mathcal{D}_s}{R} \frac {\partial \varGamma _m}{\partial \theta }\right]\left (\frac {\varGamma }{\varGamma _m}\right) \\ =2\pi R\,\sin \theta \left [u_s(\varGamma _m)\,\varGamma -\frac {\mathcal{D}_s}{R} \frac {\partial \varGamma _m}{\partial \theta }\left (\frac {\varGamma }{\varGamma _m}\right)\right]. \end{aligned} \end{equation}
Implementation of the above modifications leads to the following conservation equation for the total surface concentration of surfactant:
\begin{equation} \begin{aligned} \frac {\partial \varGamma }{\partial t} + \frac {2\varGamma \dot {R}}{R} - \frac {1}{R\sin \theta }\frac {\partial }{\partial \theta }\left [\sin \theta \, u_s(\varGamma _m)\,\varGamma -\sin \theta \,\frac {\mathcal{D}_s}{R} \frac {\partial \varGamma _m}{\partial \theta }\left (\frac {\varGamma }{\varGamma _m}\right)\right] = j_b, \end{aligned} \end{equation}
with the exchange rate with the bulk given by (2.11). It is noted that the internal exchange rate
$j_r$
, given by (2.12), results only in a redistribution of surfactant between the monolayer and the reservoir, and therefore does not enter the conservation equation for the total surface concentration. However, it is critical for the correct modelling of the evolution of the monolayer concentration,
$\varGamma _m$
, which is described by the following expression:
\begin{equation} \begin{aligned} \frac {\partial \varGamma _m}{\partial t} + \frac {2\varGamma _m\dot {R}}{R} - \frac {1}{R\sin \theta }\frac {\partial }{\partial \theta }\left [\sin \theta \,u_s(\varGamma _m)\,\varGamma _m-\sin \theta \,\frac {\mathcal{D}_s}{R} \frac {\partial \varGamma _m}{\partial \theta }\right] = j_b + j_r. \end{aligned} \end{equation}
2.5. Initial and boundary conditions
The evolution of liquid film thickness and surfactant surface concentration as described by (2.1), (2.18) and (2.19) – or (2.1) and (2.13) for the simpler case of weakly nonlinear oscillations – is supplemented by appropriate initial and boundary conditions. A convenient starting point for all the computations to be presented next is the equilibrium structure of a motionless alveolus. At equilibrium, both the interface surfactant concentration and the capillary pressure are uniform in
$\theta$
, and the interface attains a spherical shape with surface tension
$\sigma _{eq}=\sigma _{eq}(\varGamma _{eq})$
. Then, (2.1) has a trivial solution for the liquid film thickness,
$H(\theta)$
, which is a linear function of
$\cos \theta$
determined by two parameters, the mean thickness,
$\bar {H}$
, and the thickness over the rim,
$H(\theta _0)\equiv H_0$
(Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022). The computations with this initial condition indicate that in most cases a single oscillation cycle is sufficient to reach the periodic solution.
The existence of a continuous, and rather thick, liquid layer over the sharp rim (
$H_0=0.14\;\unicode{x03BC} \textrm {m}$
as estimated from the data of Bastacky et al. (Reference Bastacky, Lee, Goerke, Koushafar, Yager, Kenaga, Speed, Chen and Clements1995)) is attributed to the existence of repulsive disjoining pressure forces (Starov & Ivanov Reference Starov and Ivanov2004), which develop between the interface, the epithelium and surfactant aggregates squeezed in between. This picture is further supported by recent measurements (Xu et al. Reference Xu, Yang and Zuo2020) that estimate the combined aggregate and monolayer thickness to be of the order of
$120{-}150\; \textrm {nm}$
, i.e. in quantitative agreement with the data of Bastacky et al. (Reference Bastacky, Lee, Goerke, Koushafar, Yager, Kenaga, Speed, Chen and Clements1995). Apart from determining
$H_0$
, the disjoining pressure has no role in the dynamics, because its characteristic time scale is very large, even compared with the breathing frequency (Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Zelig & Haber Reference Zelig and Haber2002). The non-zero value of film thickness over the rim was shown by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) to play a key role in the development of tangential motion during the oscillation cycle.
Boundary conditions are applied at the apex,
$\theta =\pi /2$
, and along the rim of the alveolus at
$\theta =\theta _0$
. Symmetry is assumed at the alveolar apex, and thus
\begin{equation} \frac {\partial h}{\partial \theta }=\frac {\partial \varGamma }{\partial \theta }=\frac {\partial \varGamma _m}{\partial \theta }=0 \quad \textrm {at} \quad \theta =\pi /2. \end{equation}
It has been argued (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022) that the boundary conditions along the rim are critical for correct modelling of the dynamics, as they dictate the gradients with the interior that drive the flow and transport phenomena. In the present investigation, the boundary conditions proposed by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) are used, which amount to simple mass balances that equate water and surfactant accumulation over the rim to the respective fluxes from the interior at
$\theta _0$
.
In terms of the volumetric flow rate of water around the circumference of the rim,
$Q(\theta _0,t)$
, the mass balance of water takes the form
\begin{equation} \begin{gathered} -2\,Q(\theta _0,t) = \frac {d}{{\rm d}t}\left [2\pi R \sin \theta _0 \left (\pi \frac {(r_0+h_0)^2}{2}-\pi \frac {r_0^2}{2}\right)\right]. \end{gathered} \end{equation}
Similar mass balances for the total and monolayer surfactant, taking into account convection and diffusion at the interface and exchange with the bulk and the reservoir, lead to the following expressions:
\begin{equation} -2 Q_{\varGamma }(\theta _0,t) + 2\pi R \sin \theta _0 \left.\pi (r_0+h_0) j_b\right |_{\theta _0} = \frac {d}{{\rm d}t}\big [2\pi R\sin \theta _0\,\pi \left (r_0+h_0\right)\varGamma _0\big], \end{equation}
\begin{equation} -2 Q_{\varGamma m}(\theta _0,t) + 2\pi R \sin \theta _0 \left.\pi (r_0+h_0) (j_b+j_r)\right |_{\theta _0} = \frac {d}{{\rm d}t}\big [2\pi R\sin \theta _0\,\pi \left (r_0+h_0\right)\varGamma _{m0}\big]. \end{equation}
The multiplier two in. (2.21), (2.22) and (2.23) accounts for flux coming to the rim from both sides of the mid-alveolar wall and is appropriate for alveoli in close packing (type A alveoli), as they occur in the lowest respiratory bronchioles. The minus sign indicates flow in the negative
$\theta \hbox{-}$
direction, i.e. towards the rim.
The application of mass balances only – and not flow equations – over the rim is based on the assumption that the rim dynamics does not influence the problem, or, in other words, that it is totally enslaved to the dynamics of the alveolus. This assumption has been justified by an order-of-magnitude analysis, which revealed that Marangoni forcing is much faster than capillary drainage (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022). Thus, the flow induced inside the alveolus by Marangoni stresses dominates the system dynamics.
3. Scaling and numerical solution
3.1. The final equations in dimensionless form
The characteristic scales used to non-dimensionalise the problem variables are taken from the equilibrium conditions. Thus, we consider a motionless alveolar cap of radius
$\bar {R}$
, coated by a liquid film whose mean thickness is
$\bar {H}$
. With these choices, the lubrication parameter is formally defined as
\begin{equation} \epsilon =\frac {\bar {H}}{\bar {R}}. \end{equation}
The liquid is loaded by surfactant aggregates and equilibrates with an adsorbed monolayer of surface concentration
$\varGamma _{\textit{eq}}$
, which results in surface tension
$\sigma _{eq}$
and surface elasticity
$E_{eq}=-\varGamma _{eq} \left. ({\rm d}\sigma /{\rm d}\varGamma) \right |_{\varGamma =\varGamma _{eq}}$
. Finally, the characteristic time is the breathing period,
$T$
.
With the above scales, the following dimensionless variables are defined:
$h^*=h/\bar {H}$
,
$H^*=H/\bar {H}$
,
$R^*=R/\bar {R}$
,
$r_0^*=r_0/\bar {R}$
,
$t^*=t/T$
,
$\varGamma ^*=\varGamma /\varGamma _{eq}$
,
$\sigma ^*=\sigma /\sigma _{eq}$
,
$E^*=E/E_{eq}$
,
$p^*=p\,\bar {R}^2/(\sigma _{eq}\bar {H})$
,
$j_b^*=j_b T/\varGamma _{eq}$
and
$j_r^*=j_r T/\varGamma _{eq}$
. The governing evolution equations are transformed as follows in terms of the dimensionless variables:
\begin{equation} \frac {\partial h^*}{\partial t^*} + \frac {2h^*\dot {R^*}}{R^*} + \frac {1}{R^{*2}\sin \theta } \frac {\partial Q^*}{\partial \theta } = 0,\end{equation}
\begin{equation} \frac {\partial \varGamma ^*}{\partial t^*} + \frac {2\varGamma ^*\dot {R^*}}{R^*}+ \frac {1}{R^{*2}\sin \theta }\frac {\partial Q_{\varGamma }^*}{\partial \theta }-j_b^* = 0,\end{equation}
and
\begin{equation} \frac {\partial \varGamma _m^*}{\partial t^*} + \frac {2\varGamma _m^*\dot {R^*}}{R^*}+ \frac {1}{R^{*2}\sin \theta }\frac {\partial Q_{\varGamma m}^*}{\partial \theta }-(j_b^*+j_r^*) = 0. \end{equation}
Terms
$Q^*,\;Q_{\varGamma }^*, \;Q_{\varGamma m}^*$
are respectively the dimensionless volumetric flow rates of water, and total and monolayer surfactant, the former scaled by
$\bar {Q}=2\pi \bar {R}^2\bar {H}/T$
and the two latter by
$\bar {Q}_{\varGamma,eq}=2\pi \bar {R}^2 \varGamma _{eq}/T$
.
Three trivial algebraic manipulations are undertaken in order to derive the final evolution equations from ((3.2)–(3.4)). First, the independent spatial variable
$\theta$
is replaced by
$x=-\cos \theta$
. Second, the fluxes of water and surfactant are expanded according to the expressions in (2.1), (2.15) and (2.17), using also (2.2) and (2.16) for the pressure and the surface velocity. Third, to avoid symbol overload, the dimensionless variables are written from now on without the asterisk, e.g.
$h$
signifies
$h^*$
,
$\varGamma$
signifies
$\varGamma ^*$
,
$R$
signifies
$R^*$
etc. There should be no confusion with this last change, since only dimensionless variables appear in the rest of the paper.
Thus, the final evolution equations are
\begin{equation} \frac {\partial h}{\partial t} + \frac {2h\dot {R}}{R} + \frac {1}{R^2} \frac {\partial Q}{\partial x} = 0,\end{equation}
\begin{equation} \frac {\partial \varGamma }{\partial t} + \frac {2\varGamma \dot {R}}{R}+ \frac {1}{R^2}\frac {\partial Q_{\varGamma }}{\partial x}-j_b = 0,\end{equation}
and
\begin{equation} \frac {\partial \varGamma _m}{\partial t} + \frac {2\varGamma _m\dot {R}}{R}+ \frac {1}{R^2}\frac {\partial Q_{\varGamma m}}{\partial x}-(j_b+j_r) = 0,\end{equation}
which are solved for initial conditions
\begin{equation} h(x,0)=H(x), \quad \varGamma (x,0)=\varGamma _m(x,0)=1,\end{equation}
and subject to the following boundary conditions at
$x_0=-\cos \theta _0$
, with subscript
$0$
signifying value at the rim:
\begin{equation} -2\,Q(x_0,t)=\sqrt {1-x_0^2}\;\pi R\left [\left (r_0+\epsilon h_0\right)\frac {{\rm d}h_0}{{\rm d}t}+h_0\left (r_0+\epsilon \frac {h_0}{2}\right)\frac {\dot {R}}{R}\right],\end{equation}
\begin{equation} -2\,Q_{\varGamma }(x_0,t)=\sqrt {1-x_0^2}\;\pi R\left [(r_0+\epsilon h_0)\left (\frac {{\rm d}\varGamma _0}{{\rm d}t}-j_{b0}\right)+\epsilon \varGamma _0\frac {{\rm d}h_0}{{\rm d}t}+\varGamma _0(r_0+\epsilon h_0)\frac {\dot {R}}{R}\right],\end{equation}
\begin{align}\notag &-2\,Q_{\varGamma m}(x_0,t)\\ &\quad =\sqrt {1-x_0^2}\;\pi R\left [(r_0+\epsilon h_0)\left (\frac {{\rm d}\varGamma _{m0}}{{\rm d}t}-j_{b0}-j_{r0}\right)+\epsilon \varGamma _{m0}\frac {{\rm d}h_0}{{\rm d}t}+\varGamma _{m0}(r_0+\epsilon h_0)\frac {\dot {R}}{R}\right]. \end{align}
The symmetry boundary conditions at
$\theta =\pi$
are now trivially satisfied, provided the derivatives
$\partial h/\partial x,\;\partial \varGamma /\partial x,\;\partial \varGamma _m/\partial x$
are finite at
$x=1$
.
The fluxes that appear in ((3.5)–(3.7)) are given in expanded form by the following expressions:
\begin{equation} Q=\epsilon ^3\frac {Ca^{-1}}{3R^4}\sigma (\varGamma _m) h^3 \big(1-x^2 \big) \frac {\partial }{\partial x}(2h+g)-\epsilon \frac {Ma}{2R^2}h^2 \big(1-x^2 \big) \frac {E(\varGamma _m)}{\varGamma _m}\frac {\partial \varGamma _m}{\partial x},\end{equation}
\begin{equation} \begin{aligned} Q_{\varGamma }&=\epsilon ^3\frac {Ca^{-1}}{2R^4}\varGamma \sigma (\varGamma _m) h^2 \big(1-x^2 \big) \frac {\partial }{\partial x}(2h+g) - \epsilon \frac {Ma}{R^2}\varGamma h \big(1-x^2 \big) \frac {E(\varGamma _m)}{\varGamma _m}\frac {\partial \varGamma _m}{\partial x} \\ &\quad-\frac {Pe_s^{-1}}{R^2} \big(1-x^2 \big)\frac {\partial \varGamma _m}{\partial x}\left (\frac {\varGamma }{\varGamma _m}\right),\end{aligned} \end{equation}
\begin{equation} \begin{aligned} Q_{\varGamma m}&=\epsilon ^3\frac {Ca^{-1}}{2R^4}\varGamma _m\,\sigma (\varGamma _m) h^2 \big(1-x^2 \big) \frac {\partial }{\partial x}(2h+g) \\ &\quad -\epsilon \frac {Ma}{R^2}\varGamma _m\, h \big(1-x^2 \big) \frac {E(\varGamma _m)}{\varGamma _m}\frac {\partial \varGamma _m}{\partial x}- \frac {Pe_s^{-1}}{R^2}\big(1-x^2 \big)\frac {\partial \varGamma _m}{\partial x},\end{aligned} \end{equation}
\begin{equation} j_b = St_b \left [ (\varGamma _{m,\infty }-\varGamma _m)-\frac {\varGamma _m}{KC_{10}}\right],\end{equation}
\begin{equation} j_r = St_r (\varGamma _{m,\infty }-\varGamma _m) (\varGamma -\varGamma _m) - St_{rc}\frac {1}{\left (\varGamma _{max}-\varGamma _m\right)^2},\end{equation}
where, as a reminder, the properties
$\sigma$
and
$E$
in (3.12)–(3.14) are written with the dependence on monolayer concentration made explicit. It is noted that the function
$g(x,t)$
, which appears in these expansions, is a substitution for the expression
\begin{equation} g(x,t)=\frac {\partial }{\partial x}\left (\big(1-x^2 \big)\frac {\partial h}{\partial x}\right). \end{equation}
This substitution is necessary for the finite-element solution of the problem, as the introduction of the new unknown function
$g(x,t)$
eliminates the higher than second-order derivatives in
$h$
. Finally, the dimensionless numbers that appear in (3.12)–(3.16) are defined as follows:
\begin{equation} \begin{aligned} Ca^{-1}=\frac {\sigma _{eq}}{\mu (\bar {R}/T)}\;, \quad Ma=\frac {E_{eq}}{\mu (\bar {R}/T)}\;, \quad Pe_s^{-1}=\frac {\mathcal{D}_s T}{\bar {R}^2}\;,\\ St_b=\frac {T}{1/(k_{ads}C_{10})}\;, \quad St_r=\frac {T}{1/(k_r\varGamma _{eq})}\;, \quad St_{rc}=\frac {T}{\varGamma _{eq}^3/k_{rc}}, \end{aligned} \end{equation}
and may be viewed as ratios of the breathing period to the characteristic times of capillary, Marangoni, diffusion, adsorption and replenishment phenomena.
3.2. Numerical solution
The evolution equations ((3.5)–(3.7)) are discretised in space with a Galerkin finite-element method. The nodes of the independent spatial variable
$x \in [x_0,1]$
are clustered close to
$x_0$
, where the solution changes faster. Within each element, the primary unknowns,
$h,\,g,\,\varGamma,\,\varGamma _m$
, are expanded by quadratic Lagrangian basis functions,
$\phi _i$
, and the evolution equations are weighted integrally with the same basis functions. The final weak forms are as follows, after applying integration by parts to reduce the order of differentiation:
\begin{equation} \int _{x_0}^{1} \frac {\partial h}{\partial t}\phi _i \,{\rm d}x + \int _{x_0}^{1} \frac {2h\dot {R}}{R}\phi _i \,{\rm d}x + \frac {1}{R^2}\left [Q\phi _i\right]_{x_0}^{1} -\frac {1}{R^2}\int _{x_0}^{1} Q\frac {d\phi _i}{{\rm d}x}\,{\rm d}x = 0,\end{equation}
\begin{equation} \int _{x_0}^{1} g\,\phi _i \,{\rm d}x - \left [\big(1-x^2 \big)\frac {\partial h}{\partial x}\phi _i\right]_{x_0}^{1} + \int _{x_0}^{1} \big(1-x^2 \big)\frac {\partial h}{\partial x}\frac {d\phi _i}{{\rm d}x}\,{\rm d}x = 0,\end{equation}
\begin{equation} \int _{x_0}^{1} \frac {\partial \varGamma }{\partial t}\phi _i \,{\rm d}x + \int _{x_0}^{1} \frac {2\varGamma \dot {R}}{R}\phi _i \,{\rm d}x + \frac {1}{R^2}\left [Q_{\varGamma }\phi _i\right]_{x_0}^{1} -\frac {1}{R^2}\int _{x_0}^{1} Q_{\varGamma }\frac {d\phi _i}{{\rm d}x}\,{\rm d}x-\int _{x_0}^{1}j_b\phi _i \,{\rm d}x= 0,\end{equation}
and
\begin{equation} \begin{aligned} \int _{x_0}^{1} \frac {\partial \varGamma _m}{\partial t}\phi _i \,{\rm d}x + \int _{x_0}^{1} \frac {2\varGamma _m\dot {R}}{R}\phi _i \,{\rm d}x + \frac {1}{R^2}\left [Q_{\varGamma m}\phi _i\right]_{x_0}^{1} -\frac {1}{R^2}\int _{x_0}^{1} Q_{\varGamma m}\frac {d\phi _i}{{\rm d}x}\,{\rm d}x \\ -\int _{x_0}^{1}(j_b+j_r)\phi _i \,{\rm d}x = 0.\end{aligned} \end{equation}
The fluxes in the integrated terms in (3.19)–(3.22) are by symmetry equal to zero at
$x=1$
, and are evaluated at
$x=x_0$
from the boundary conditions, (3.9)–(3.11). Time derivatives are treated with an implicit Euler scheme and the spatial integrals are evaluated by three-point Gaussian quadrature. Finally, the resulting system of nonlinear algebraic equations is solved with a Newton–Raphson iterative scheme.
4. Results
The results presented in the rest of the paper are computed using the following parameter values. The equilibrium alveolar radius and the septum thickness at the rim are
$\bar {R}=100\,\unicode{x03BC} {\rm m}$
and
$2r_0=4\,\unicode{x03BC}{\rm m}$
, respectively (Ochs et al. Reference Ochs, Nyengaard, Jung, Knudsen, Voigt, Wahlers, Richter and Gundersen2004). Alveolar oscillations are purely sinusoidal, characterised by the dimensionless amplitude,
$a$
, i.e. in dimensionless terms
\begin{equation} R(t)=1+a\sin \left (2\pi t\right).\end{equation}
The liquid layer has the density and viscosity of pure water at the physiological temperature of
$37\,^\circ\rm C$
. At equilibrium, its mean thickness is
$\bar {H}=0.3\,\unicode{x03BC}\rm m$
and its thickness over the rim is
$H_0=0.14\,\unicode{x03BC}\rm m$
(Bastacky et al. Reference Bastacky, Lee, Goerke, Koushafar, Yager, Kenaga, Speed, Chen and Clements1995; Xu et al. Reference Xu, Yang and Zuo2020). The equilibrium surface tension is taken as
$\sigma _{eq}=0.023\,{\rm N\,m}^{-1}$
, the molar surface area at zero surface pressure
$\varOmega _0=3.0\times 10^5\,\rm m^{2}\,mol^{-1}$
and the intrinsic compressibility
$\alpha =0.0052\,\rm m\,mN^{-1}$
(Wüstneck et al. Reference Wüstneck, Perez-Gil, Wüstneck, Cruz, Fainerman and Pison2005; Zuo et al. Reference Zuo, Rimei, Xianju, Jinlong, Policova and Neumann2016). Surfactant exchange between the monolayer and the bulk is quantified by the adsorption parameter
$k_{ads}C_{10}$
. Its value is zero for the case of an insoluble surfactant and takes typical values in the range
$k_{ads}C_{10}=10{-}40\,\rm s^{-1}$
for soluble surfactants (Zuo et al. Reference Zuo, Gitiafroz, Acosta, Policova, Cox, Hair and Neumann2005; Saad et al. Reference Saad, Neumann and Acosta2010; Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2021). The rate of monolayer replenishment from the reservoir during alveolar expansion is taken as
$k_r=150\,\rm m^2\,mol^{-1}\,s^{-1}$
. The values of
$\sigma _{eq}$
,
$k_{ads}C_{10},\,\alpha$
and
$k_r$
are derived from best fit with available experimental data, as described in Appendix A. The parameter
$k_{rc}$
that quantifies how abruptly the monolayer concentration arrests at the value
$\varGamma _{max}$
during monolayer collapse is set to
$k_{rc}=10^{-5}\,\rm mol^3\,m^{-6}\,s^{-1}$
.
It is noted that there is a high degree of confidence in the accuracy of the selected parameter values, for two reasons. First, each parameter affects the data (the surface tension–surface area loop) in a different way. For example
$\sigma _{eq}$
moves the entire loop up or down,
$k_{ads}C_{10}$
varies the width of the loop,
$\alpha$
modifies the slope of the compression leg and
$k_r$
restricts the upper limit of the loop. Thus, there is a single combination of best-fit parameter values. Second, the values selected ‘blindly’ by the above procedure are in very satisfactory agreement with the literature. For example, all measured data of equilibrium surface tension are in the range
$\sigma _{eq}=22{-}25\,\rm mN\,m^{-1}$
, independent from the total surfactant concentration (Wüstneck et al. Reference Wüstneck, Perez-Gil, Wüstneck, Cruz, Fainerman and Pison2005; Parra & Perez-Gil Reference Parra and Perez-Gil2015). Also, the best-fit value of the monolayer compressibility,
$\alpha$
, is in the expected range for the tilted condensed regime and – most important – in agreement with the compressibility of a pure DPPC monolayer (Schürch et al. Reference Schürch, Green and Bachofen1998; Possmayer et al. Reference Possmayer, Nag, Rodriguez, Qanbar and Schürch2001; Zhang et al. Reference Zhang, Wang, Fan and Zuo2011).
4.1. The dynamics of large-amplitude oscillations
The onset of monolayer collapse when
$\varGamma \rightarrow \varGamma _{max}$
occurs at oscillation amplitudes
$a\approx 0.05{-}0.08$
, with the exact value depending on the minimum surface tension achievable by compression,
$\sigma _{min}=\sigma (\varGamma _{max})$
, and on the surfactant solubility. The value of minimum surface tension used throughout the paper is
$\sigma _{min}=0.005\,{\rm N\,m^{-1}}$
and the dynamics of the system is first examined for an insoluble surfactant.
An overview of the behaviour of the alveolus is provided by figure 3, which shows the temporal variation of
$\varGamma$
and
$\varGamma _m$
for two oscillation cycles starting from rest at equilibrium. It is confirmed that the variation of surface surfactant concentration is mainly driven by the change in the alveolar surface area and thus is
$180^\circ$
out of phase with the oscillation amplitude. The spatial variation of
$\varGamma$
and
$\varGamma _m$
(the latter responsible for the generation of Marangoni stresses) is too small to be visually distinguished in such a plot, so the depicted behaviour is valid for all points along the alveolar interface.
The temporal variation of alveolar radius
$R(t)$
(black dashed line), monolayer concentration
$\varGamma _m$
(red dashed line) and total surface concentration
$\varGamma$
(green solid line) for oscillations with amplitude (a)
$a=0.05$
and (b)
$a=0.15$
.
Figure 3(a) refers to oscillations with an amplitude (
$a=0.05$
) which is too small to cause monolayer collapse. Thus, the curves for
$\varGamma$
and
$\varGamma _m$
are superimposed. The situation is different in figure 3(b), where the oscillation amplitude is
$a=0.15$
. During compression, the total surface concentration increases continuously while the monolayer concentration is saturated at
$\varGamma _m =1.1465$
, i.e.
$\approx 15\,\%$
higher than the equilibrium value. During re-expansion, the monolayer is enriched by transport from the reservoir, which is eventually consumed leaving once again a pure monolayer (
$\varGamma =\varGamma _m$
).
For small-amplitude oscillations, Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) developed expansions up to second order in the amplitude of the following general form for variables of interest such as surface deformation,
$h(x,t)$
, surface velocity,
$u_s(x,t)$
, volumetric flow rate,
$Q(x,t)$
and wall shear stress,
$\tau _w(x,t)$
:
\begin{equation} \varPhi (x,t)=a\,\textrm{Re}\big [\varPhi _1(x)\,e^{{i2\pi t}}\big]+ a^2\Big (\textrm{Re}\big [\varPhi _2(x)\,e^{{i4\pi t}}\big]+\varPhi _{2s}(x)\Big). \end{equation}
The analytic expressions for the amplitudes in terms of the system parameters are contained in the original publication. Given that the linear/weakly nonlinear solution is relatively simple, and that it has already been used to provide predictions beyond its formal range of validity (Livanos et al. Reference Livanos, Bouchoris, Aslani, Gourgoulianis and Bontozoglou2024), it is of interest to compare it with the full numerical solution.
Numerical results for the surface velocity are computed from (2.16), for the volumetric flow rate from (2.3) and for the wall shear stress from the following expression (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022):
\begin{equation} \tau _w=\frac {h}{R}\frac {\partial p}{\partial \theta }-\frac {1}{R}\frac {\partial \sigma }{\partial \theta }. \end{equation}
Indicative numerical results for oscillation amplitude
$a=0.0325$
are Fourier transformed and the magnitudes of the first two harmonics are compared with the weakly nonlinear predictions in figure 4(a,b). It is observed that numerics and asymptotics are in close agreement. The gross characteristics of
$h(x,t)$
,
$u_s(x,t)$
,
$Q(x,t)$
and
$\tau _w(x,t)$
as the oscillation amplitude increases, and the extent to which they may be roughly estimated from the linear prediction, are examined in more detail in Appendix B.
The first and second harmonics of the numerical solution (continuous lines) compared with the weakly nonlinear predictions (open circles).
A key characteristic of the dynamics of the present problem – already discussed in the linear limit by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) – is the contrast between the dominance of Marangoni stresses along the interface and the significant role of capillary and viscous stresses in the interior of the film. Thus, although
$Ma\gt \gt Ca^{-1}$
, all contributions must be retained for correct prediction of the flow.
The dominance of Marangoni over capillary stresses on the interface results in the instantaneous establishment of an inverse Marangoni flow (Manikantan & Squires Reference Manikantan and Squires2020) that negates the effect of variations in interfacial curvature. This behaviour is witnessed by the spatial distribution of the monolayer concentration of surfactant, which is depicted in figure 5(a,b). In order to include results for a number of representative time instants, we actually plot the difference
$(\varGamma _m(x,t)-\varGamma _m(1,t))$
, normalised with the oscillation amplitude
$a$
, so that all the curves coincide at the apex. The oscillatory shape of these curves occurs in response to variations of interfacial curvature, and is contrasted with the monotonic change of
$u_s(x,t)$
with
$x$
in figure 4(a).
The term
$[\varGamma _m(x,t)-\varGamma _m(1,t)]/a$
for time instants
$t/T=0.125$
(red),
$0.25$
(blue),
$0.375$
(green) and
$0.50$
(black) in panel (a) and
$t/T=0.625$
(red),
$0.75$
(blue),
$0.875$
(green) and
$1.0$
(black) in panel (b). The alveolar oscillations have amplitude
$a=0.05$
(dashed lines) and
$0.20$
(continuous lines).
The results in figure 5 are qualitatively similar for the relatively small oscillation amplitude
$a=0.05$
(dashed lines) and for
$a=0.20$
(continuous lines), despite the extensive monolayer collapse during compression that occurs in the second case. The behaviour indicates that the limitation imposed on
$\varGamma _m$
due to the collapse does not obstruct the spatial variations that give rise to Marangoni stresses. This happens because the necessary gradients are very small and the modelling of collapse by (2.12) allows for them to develop. Restriction of the spatial gradients of
$\varGamma _m$
and of the ensuing Marangoni stresses occurs only for extremely small values of the parameter
$k_{rc}$
, which determines the stiffness of the monolayer. This behaviour is depicted in figure 6, which shows the temporal maxima of the surface velocity,
$u_s$
, along the liquid film for different values of
$k_{rc}$
. It is observed that the Marangoni-induced surface flow is significantly affected only below
$k_{rc}=10^{-6}$
. The results presented in the rest of the paper are computed using
$k_{rc}=10^{-5}$
and the role of
$k_{rc}$
is not further investigated, given the apparent lack of measured data.
The maxima of surface velocity for oscillations with amplitude
$a=0.20$
and monolayer stiffness parameter
$k_{rc}=10^{-4},\,10^{-5},\,10^{-6},\,10^{-7}$
and
$10^{-8}$
(red, black dashed, green, blue dashed and purple lines, respectively).
Although the main flow variables appear to scale roughly with the oscillation amplitude, there are important modifications with increasing nonlinearity and in particular with the onset of monolayer collapse. A first impression is provided by considering the tangential velocity field inside the liquid layer. Figure 7(a,b) shows
$(u_{\theta }/a)$
for equidistant time instants during a breathing cycle for
$a=0.01$
and
$a=0.20$
, respectively. The similarity of figure 7(a) to an inverted Stokes layer is evident, with the main difference being that the motion is now triggered not by tangential oscillation of the wall but by the effect of Marangoni stresses. Another important difference is that momentum transfer from the surface towards the wall is affected not only by viscosity but also by capillary pressure, induced on the free surface by local variation of curvature.
The profile of tangential velocity,
$u_{\theta }$
, normalised with the oscillation amplitude
$a$
, inside the liquid layer, plotted for different time instants during a breathing cycle: (a)
$a=0.01$
and (b)
$a=0.20$
.
Recalling that Marangoni stresses are triggered by the periodic variation of the monolayer concentration inside the alveolus, the symmetric shape of figure 7(a) is explained by the fact that, in the linear limit, the variation in wall surface area is proportional to the oscillation amplitude. Thus, the change in surface tension is symmetric during expansion and contraction of the alveolus. However, this symmetry breaks with increasing oscillation amplitude as
$O(a^2)$
effects come into play. Then, the decrease in surface area during contraction is less than the increase during expansion, and capillary pressure introduces asymmetry. The effect intensifies much more with the onset of monolayer collapse, because then the surface tension, instead of varying in response to the change in surface area, remains constant at
$\sigma _{min}=\sigma (\varGamma _{max})$
. The asymmetric variation of the velocity profile for
$a=0.20$
is depicted in figure 7(b).
Sensitive indicators of the effects of viscosity and capillarity on the transmission of momentum through the liquid film are the volumetric flow rate,
$Q$
, and the wall shear stress,
$\tau _w$
. Thus, figures 8(a,b) and 9(a,b) respectively plot data for the temporal variation of
$Q$
and
$\tau _w$
, during one breathing cycle. Figures 8(a) (for
$Q$
) and 9(a) (for
$\tau _w$
) refer to a single location close to the alveolar rim at
$\theta \approx 37^\circ$
. Black lines correspond to amplitudes
$a=0.01{-}0.065$
that do not reach monolayer collapse, and red lines to amplitudes
$a=0.08{-}0.20$
that cause collapse. In the linear limit (continuous black lines) both variations of
$Q$
and
$\tau _w$
are sinusoidal, i.e. perfectly symmetric in time. Asymmetry develops gradually with increasing amplitude in the compression part of the cycle and, with the onset of collapse, the flow rate abruptly decelerates and then for a short time changes direction. Similar abrupt change (a spike of short duration) is observed in figure 9(a) in the wall shear stress. The disturbances in
$Q$
and
$\tau _w$
intensify with the increase in oscillation amplitude, but magnification of the data shows that there is no discontinuity in slope.
Figures 8(b) and 9(b) depict respectively the temporal variations of
$Q$
and
$\tau _w$
for a number of locations
$\theta \approx 33^\circ,\,37^\circ,\,42^\circ,\,56^\circ$
and
$74^\circ$
and show that the spikes are more intense close to the rim, but are also observed, with diminishing intensity, deeper inside the alveolus. The strong effect of nonlinearity may be appreciated by recalling that, in the linear limit, all the curves in figures 8(b) and 9(b) would be simple sinusoids of different amplitudes.
(a) The temporal variation of volumetric flow rate at
$\theta \approx 37^\circ$
during one breathing cycle for
$k_{ads}C_{10}=0\,{\rm s}^{-1}$
. Black lines are for
$a=0.01$
(continuous),
$0.0325$
(dashed),
$0.05$
(dash-dotted) and
$0.065$
(dotted). Red lines are for
$a=0.08$
(continuous),
$0.10$
(dashed),
$0.15$
(dash-dotted) and
$0.20$
(dotted). (b) The temporal variation of volumetric flow rate for
$a=0.20$
at locations
$\theta \approx 33^\circ$
(dashed),
$37^\circ$
(continuous),
$42^\circ$
(dash-dotted),
$56^\circ$
(dotted) and
$74^\circ$
(thick dots).
(a) The temporal variation of wall shear stress at
$\theta \approx 37^\circ$
during one breathing cycle for
$k_{ads}C_{10}=0\,{\rm s}^{-1}$
. Black lines are for
$a=0.01$
(continuous),
$0.0325$
(dashed),
$0.05$
(dash-dotted) and
$0.065$
(dotted). Red lines are for
$a=0.08$
(continuous),
$0.10$
(dashed),
$0.15$
(dash-dotted) and
$0.20$
(dotted). (b) The temporal variation of wall shear stress for
$a=0.20$
at locations
$\theta \approx 33^\circ$
(dashed),
$37^\circ$
(continuous),
$42^\circ$
(dash-dotted),
$56^\circ$
(dotted) and
$74^\circ$
(thick dots).
The sudden variation in wall shear stress during deep breathing, which is related to the onset of monolayer collapse and has a short time scale compared with the breathing period, appears to increase significantly the complexity in the temporal variation of shear forces on the epithelium compared with a simple sinusoid. In addition, the expected random variation in oscillation amplitude during consecutive breathing cycles will further increase the shear stress variability both in magnitude and in phase. This behaviour may have an interesting physiological effect, as discussed in the Concluding Remarks.
The analysis up to this point referred to insoluble surfactant. For soluble surfactant with
$k_{ads}C_{10}\sim 10{-}40\,\rm {\rm s}^{-1}$
, the characteristic time of exchange between the monolayer and the bulk is
$t_{ads}\sim 10^{-1}{-}10^{-2}\, s$
, which is much slower than the characteristic times of capillary and Marangoni flows, cf. (3.18). As a result, solubility has negligible effect at leading order and the results for the main flow parameters presented above retain the same qualitative character with relatively minor quantitative changes. Thus, they are not deemed worth reproducing. However, surfactant solubility has a significant effect on second-order phenomena to be considered next, namely the slow displacement of nanoparticles deposited on the liquid layer and the drift of surfactant through the alveolar opening.
4.2. Alveolar clearance of insoluble particles
Small particles contained in the inhaled air are transported distally into the lungs and many of them deposit. Clearance of insoluble particles, deposited from the trachea down to the respiratory bronchioles, occurs relatively fast (a large fraction is cleared within a few hours) through the well-studied mucociliary escalator mechanism (Smith, Eamonn & Blake Reference Smith, Eamonn and Blake2008). Alveolar clearance is many orders of magnitude slower and its mechanism is still open to investigation. Experiments with human volunteers who inhaled
${\sim} 1\,\unicode{x03BC}\rm m$
diameter fused aluminosilicate particles (Bailey et al. Reference Bailey, Frances and James1982, Reference Bailey, Fry and James1985) showed that the particles, which were retained after six days, cleared at a rate described by a two-component exponential function, the first component with a half-time of tens of days and the second with a half-time of hundreds of days.
Potential routes of alveolar clearance are surface transport from the alveoli to the ciliated regions of the lung and/or translocation through the epithelium and secretion by the lymphatic system (Lippmann, Yeates & Albert Reference Lippmann, Yeates and Albert1980; Oberdoerster Reference Oberdoerster1988). Phagocytosis by alveolar macrophages is a major mechanism of large particle arrest, although the macrophages must also be eventually transported to ciliated regions. However, phagocytosis appears to be far less efficient for particles of nanometre size (Oberdoerster Reference Oberdoerster1988; Ferin Reference Ferin1994). The role of surface tension, in combination with its variations induced by the breathing manoeuvres, has been hypothesised long ago as the main transport mechanism (Kilburn Reference Kilburn1968). However, this hypothesis could not be tested further because of the limited understanding of the fluid dynamics involved.
The present model provides an opportunity to investigate whether the flow inside the alveolus has a clearance effect. To this end, light, insoluble nanoparticles are assumed to stay on the surface of the liquid layer and to behave as passive tracers. Their net displacement per breathing cycle is thus computed by integrating in time the equation of motion in the lubrication limit,
$u_{\theta,s}\,{\rm d}t=R\,d\theta$
, by a fourth-order Runge–Kutta scheme. The results are expressed as displacement angle from the initial position, normalised with the square of the amplitude, and are presented in figure 10(a,b) for insoluble surfactant (continuous lines) and for surfactants with
$k_{ads}C_{10}=12\,\rm {\rm s}^{-1}$
(dashed lines) and
$k_{ads}C_{10}=28\,\rm {\rm s}^{-1}$
(dashed–dotted lines).
The displacement of a floating nanoparticle after one breathing cycle, as function of the initial deposition location, computed for surfactant solubilities
$k_{ads}C_{10}=0\,{\rm s}^{-1}$
(continuous lines),
$k_{ads}C_{10}=12\,{\rm s}^{-1}$
(dashed lines) and
$k_{ads}C_{10}=28\,{\rm s}^{-1}$
(dashed–dotted lines). (a) Small oscillation amplitudes
$a=0.0325$
(red),
$0.05$
(blue) and
$0.065$
(green). (b) Large oscillation amplitudes
$a=0.10$
(red),
$0.15$
(blue) and
$0.20$
(green).
Figure 10(a) refers to finite but relatively small oscillation amplitudes, i.e. calm breathing. The variation with amplitude appears as an
$O(a^2)$
effect, with this dependence slightly weakening with increasing
$a$
. The most interesting findings are that (i) under all conditions, particle drift is negative, i.e. all particles move toward the rim, and (ii) that increased surfactant solubility strongly intensifies the clearance rate. In general, the drift is more intense in the first
$50^\circ$
from the alveolar opening and maximises right on the rim. The local maximum very close to the rim, which is more pronounced at high solubilities and corresponds to a minimum in the absolute value of drift, is expected to contribute to the local accumulation of particles.
Figure 10(b) refers to large-amplitude oscillations, i.e. deep breathing, and the strong positive effect of surfactant solubility on clearance rate is again confirmed. However, there are two major differences from figure 10(b). First, there is now a positive effect of oscillation amplitude (beyond the
$O(a^2)$
dependence), which intensifies with solubility. Second, with increasing oscillation amplitude the displacement at the rim decreases and in some cases (small solubility, large amplitude) becomes even positive. This means that, in those cases, particles at the rim may be temporarily obstructed by the flow from entering the alveolar ducts. However, given that particle displacements are still very small and thus clearance time scales remain of the order of many tens or hundreds of days, the high-amplitude results are not considered any further (no one can breath deeply all day long and for many days).
The information contained in figure 10(a) may be utilised directly to compute the locations of passive tracers after a large number of breathing cycles. Figure 11(a) shows the drift of particles, initially equi-distributed azimuthally, during the passage of
$300$
days, when
$k_{ads}C_{10}=28$
and
$a=0.0325$
(
$\pm 10\,\%$
variation of acinar volume from the rest value). The fraction of particles that have reached the rim (considered cleared from the alveolus) is plotted in figure 11(b) versus time (circles) and is shown to obey an exponential variation (red line) with half-life
$\tau _{1/2}=301$
days. This value appears very small when compared with the time scale of 1–2 days for clearance of the conducting airways (nasal passages and tracheobronchial tree). However, it is in quantitative agreement with the results of in vivo experiments with insoluble particles by Bailey et al., which gave alveolar clearance half-times of 320 (Bailey, Frances & James Reference Bailey, Frances and James1982) and 350 days (Bailey, Fry & James Reference Bailey, Fry and James1985).
(a) The trajectory with time of particles initially equi-distributed inside the alveolus, for surfactant solubility
$k_{ads}C_{10}=28\,{\rm s}^{-1}$
and oscillation amplitude
$a=0.0325$
. (b) The fraction of particles that have reached the rim versus time (blue circles) and an exponential fit with half-life 301 days (red line).
A second exponential decay with much shorter half-time (
$\tau _{1/2}\sim 10^1$
) is expected to occur if, on top of the uniform distribution, there is an additional amount of particles deposited close to the alveolar entrance. Indeed, there is experimental evidence for such behaviour (Zeltner et al. Reference Zeltner, Sweeney, Skornik, Feldman and Brain1991), which was also mechanistically interpreted as a result of interaction of the ductal airflow with recirculating flow inside the alveolus (Henry & Tsuda Reference Henry and Tsuda2010; Tsuda, Henry & Butler Reference Tsuda, Henry and Butler2013). This tentative mechanism is valid only for the proximal alveolar generations, as it is believed that recirculation is absent in the most distal generations and air enters the alveolus with roughly radial direction (Kumar et al. Reference Kumar, Tawhai, Hoffman and Lin2009; Darquenne, Harrington & Prisk Reference Darquenne, Harrington and Prisk2009).
4.3. Surfactant drift through the alveolar opening
The weakly nonlinear analysis of Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) showed that, for soluble surfactants, surfactant flux
$Q_{\varGamma }$
has a steady streaming component, which results in an exit flux,
$Q_{\varGamma 0}$
, from the alveolar opening. This prediction is of particular interest in relation to techniques that probe exhaled air for micro-droplets originating from distal lung areas (Shmyrov et al. Reference Shmyrov, Mizev, Mizeva and Shmyrova2021). The accepted view is that these droplets form during re-opening of the small airways, after their closure during exhalation down to residual volume (Almstrand et al. Reference Almstrand, Bake, Ljungström, Larsson, Bredberg, Mirgorodskaya and Olin2010; Grotberg Reference Grotberg2011). Airway lining fluid is known to contain surfactant molecules, whose analysis indicates that they most probably originate from the alveoli, where they are actually produced (Bernhard et al. Reference Bernhard, Haagsman, Tschernig, Poets, Postle, van Eijk and von der Hardt1997). Thus, the present prediction of a steady drift from the alveoli represents a potential mechanism for the transport of surfactant to the airways.
The dependence of the steady drift,
$Q_{\varGamma 0}$
, on surfactant solubility,
$k_{ads}C_{10}$
, is investigated in detail for a wide range of oscillation amplitudes in figure 12(a,b). The weakly nonlinear prediction
$Q_{\varGamma 0}/a^2$
is plotted in figure 12(a) as a continuous red line and the drift is observed to increase from zero for the insoluble surfactant to a maximum of
$Q_{\varGamma 0}\approx 0.0574\,a^2$
at
$k_{ads}C_{10}=22\,{\rm s}^{-1}$
, and then slowly decrease for higher solubilities. The results for oscillations of finite amplitude indicate a gradual decrease in the coefficient of the
$O(a^2)$
variation and a shift of the maximum to higher values. Thus, figure 12(b) (which is a mere re-plot of the data in figure 12
a) shows that
$Q_{\varGamma 0}$
increases roughly linearly with
$a$
. For the same amplitude, the increase with solubility is a consequence of the fact that the data are below the respective maximum for each amplitude, as depicted in figure 12(a).
(a) The dependence of dimensionless surfactant drift on solubility, the former normalised with the square of the oscillation amplitude. Weakly nonlinear prediction (red line) and results for
$a=0.05$
(blue),
$0.10$
(green),
$0.15$
(red) and
$0.20$
(black circles). (b) The variation of dimensionless surfactant drift,
$Q_{\varGamma 0}$
, with amplitude for solubility
$k_{ads}C_{10}=12\,{\rm s}^{-1}$
(red crosses),
$28\,{\rm s}^{-1}$
(blue squares) and
$40\,{\rm s}^{-1}$
(green squares).
5. Concluding remarks
An oscillating spherical cap, lined internally with a surfactant-laden liquid film, is considered a model of the dynamics of a single alveolus. The linear and weakly nonlinear analysis of this system around equilibrium by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) – in the lubrication limit of quasi-steady Stokes flow – is presently extended numerically to large-amplitude oscillations. The same boundary conditions are used at the rim, which enforce local conservation of the mass of water and surfactant without imposing any artificial constraint on the magnitude of the film thickness and surfactant concentration. In the case of a soluble surfactant, the adsorption to the liquid–air interface is described by Langmuir kinetics, modified by the inclusion of the intrinsic compressibility of the adsorbed monolayer. The parameters of the model of surfactant dynamics are best-fitted to achieve quantitative agreement with the behaviour of actual pulmonary surfactant preparations, as described in Appendix A.
The main qualitative difference from previous analysis, caused by large-amplitude oscillations, is the inclusion of phenomena related to monolayer collapse and replenishment. It is known (Krueger & Gaver Reference Krueger and Gaver2000; Lee Reference Lee2008; Parra & Perez-Gil Reference Parra and Perez-Gil2015) that during strong compression of the interface, the surfactant monolayer reaches a limiting surface concentration,
$\varGamma _{max}$
, corresponding to minimum surface tension,
$\sigma _{min}=\sigma (\varGamma _{max})$
, and then collapses to form a sub-surface reservoir consisting of protrusions toward the interior that remain linked to the monolayer. This reservoir replenishes the monolayer during re-expansion of the interface and thus restricts the rise in surface tension.
The description of monolayer collapse and replenishment necessitates differentiating between the monolayer concentration,
$\varGamma _m$
, and the total surface concentration,
$\varGamma$
(with
$\varGamma -\varGamma _m$
being the reservoir concentration), and independent conservation equations are now considered for
$\varGamma _m$
and for
$\varGamma$
. A key addition is the introduction of a novel kinetic expression to describe the ‘two-way’ exchange of surfactant between the monolayer and the reservoir. The kinetic expression contains a parameter that quantifies the stiffness of the monolayer, i.e. the approach to the maximum concentration needed to initiate collapse.
It is confirmed that the flow is driven by Marangoni stresses, developed as a result of differences in surfactant concentration between the rim and the interior of the alveolus. It is also found that the collapse of the monolayer does not hinder significantly the development of these stresses, unless the stiffness parameter becomes near singular. Thus, the first-order behaviour is similar for all oscillation amplitudes and a rough estimate of magnitude for the main flow variables can be provided from linear theory (see Appendix B below).
Capillary stresses, although dominated by Marangoni stresses on the interface, have a significant effect on internal variables such as the volumetric flow rate through the film and the shear stress along the epithelial wall. As capillary stresses are strongly modified during collapse because surface tension remains constant with further compression, the temporal variation of
$Q$
and
$\tau _w$
deviates drastically from linear behaviour. Important second-order effects, such as the slow displacement of passive tracers on the interface and the drift of surfactant from the alveolus are investigated in detail and found to depend on both oscillation amplitude and surfactant solubility.
It is argued that some of the aforementioned results may tentatively be of physiological significance. The predicted nonlinear behaviour induces abrupt changes in wall shear stress, which vary both in magnitude and phase with oscillation amplitude. The overall effect, including the expected random variation of oscillation amplitude during consecutive breathing cycles, is equivalent to the imposition of near-random variability in the epithelial shear. As was demonstrated in vitro by Arold, Bartola & Suki (Reference Arold, Bartolák-Suki and Suki2009), stretching with random temporal variability, exerted on type-II alveolar cells, leads to enhanced surfactant secretion in comparison with purely periodic stretching. However, the rate of surfactant drift from the alveolus is presently predicted to increase with oscillation amplitude. Thus, the combination of the two effects appears to contribute to homeostasis at the alveolar level.
Also, the present prediction for the displacement of passive tracers by Marangoni stresses indicates that nano-particles deposited on the liquid layer may be cleared from the alveolus by the oscillation-induced flow field. The predicted exit half-time appears to agree with in vivo data, lending support to the argument that Marangoni flow provides a significant mechanism of alveolar clearance. However, the data did not involve direct measurement of alveolar retention, so other possibilities, such as uptake by alveolar macrophages and excretion involving the inflammatory system and lymphatics, cannot be excluded.
However, the passive tracer assumption is gradually invalidated with increasing particle diameter and two-way coupling between the particles and the flow will start to influence the dynamics. Such a dynamics would be more accurately described by interface-resolved simulations (Uhlmann Reference Uhlmann2008; Costa, Brandt & Picano Reference Costa, Brandt and Picano2020), which, however, are orders of magnitude more demanding in terms of computer resources and would also require a radical reformulation of the problem. Interface-resolved simulations would also allow examination of the influence of the surfactant aggregates on the dynamics of the liquid film, an effect that has been currently neglected.
The present investigation is based on the assumed alveolar shape of a spherical cap in self-similar oscillation. Although this is the most frequently used model in the literature, both the true shape and the mode of oscillation of alveoli are subject to extensive discussion (Smaldone & Mitzner Reference Smaldone and Mitzner2012; Sznitman Reference Sznitman2022). On the one hand, space-filling constraints of the acinus promote polyhedral over spherical shapes. On the other hand, the self-similar oscillation, which corresponds to isotropic expansion i.e. alveolar surface area proportional to alveolar volume raised to the
$2/3$
power, is but one of the possible modes of oscillation. There are other possibilities (Smaldone & Mitzner Reference Smaldone and Mitzner2012), with most exotic one being alveolar recruitment, at least for part of the acinus. Recruitment assumes an accordion-like behaviour, with the alveolus collapsed during deep exhalation and reforming during inhalation. Thus, the alveolar volume changes while the epithelial surface area remains roughly constant.
Despite the above uncertainty in alveolar shape and oscillation pattern, it is believed that the conclusions of the present study are of more general validity, at least qualitatively. In particular, it is argued that periodic changes in interfacial surface area occur during the breathing manoeuvre for all the tentative expansion/contraction patterns, and a schematic demonstration for the ‘accordion’ model is shown in figure 13. Thus, gradients in surfactant concentration between the rim and the interior of the alveolus, which constitute the focus of the present study, remain in all cases the driving force for the shearing flow of the liquid layer that lines the alveolar epithelium.
A schematic indication of the change in interfacial surface area during recruiting of a contracted alveolus.
Declaration of interests
The author reports no conflict of interest.
Appendix A.
The model of surfactant dynamics described in § 2.3 has been validated with data from Saad et al. (Reference Saad, Neumann and Acosta2010) and Xu et al. (Reference Xu, Yang and Zuo2020). These authors measured the dynamic surface tension of a drop of pulmonary surfactant preparation, which experiences a periodic change in its surface area
$A$
. The drop rests on a flat circular substrate with sharp edges, so that its surface remains spherical, and the variations in surface area are caused by controlled rates of inflow and outflow of liquid from a hole at the centre of the substrate.
Under these conditions, surfactant concentrations,
$\varGamma$
and
$\varGamma _m$
, vary only with time, and the respective conservation equations simplify to the following expressions in terms of the known variables
$A$
and
${\rm d}A/{\rm d}t$
:
\begin{align} \frac {{\rm d}\varGamma }{{\rm d}t} + \frac {\varGamma }{A}\frac {{\rm d}A}{{\rm d}t} + k_{ads}C_{10}\left [(\varGamma _{m,\infty }-\varGamma _m)-\frac {\varGamma _m}{K\,C_{10}} \right] = 0, \\[-29pt] \nonumber \end{align}
\begin{align} &\frac {{\rm d}\varGamma _m}{{\rm d}t} + \frac {\varGamma _m}{A}\frac {{\rm d}A}{{\rm d}t} + k_{ads}C_{10}\left [(\varGamma _{m,\infty }-\varGamma _m)-\frac {\varGamma _m}{K\,C_{10}} \right] \nonumber\\ &\quad - k_r \left ( \varGamma _{m,\infty } - \varGamma _{m} \right)\left (\varGamma - \varGamma _m\right) + \frac {k_{rc}}{\left (\varGamma _{max}-\varGamma _m\right)^2} = 0. \end{align}
Equations (A1)–(A2) are integrated with a Newton–Raphson iterative scheme, using implicit Euler differencing for the derivatives.
The key parameters of the surfactant,
$KC_{10},\,k_{ads}C_{10},\,\alpha,\,k_r$
and
$k_{rc}$
, are fitted, so that predictions by the above model agree quantitatively with the data. Comparisons are shown in figure 14(a–d) and the best-fit values are physically meaningful and in the expected range (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2021). It is noted that some of these values differ from the ones reported by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2021), because an improved model is presently used to describe monolayer collapse during deep compression and replenishment during re-expansion.
Predictions (red lines) and data (blue circles) by Saad et al. (Reference Saad, Neumann and Acosta2010) of the temporal variation of surface tension during periodic compression and expansion of the free surface of a droplet doped with surfactant. Panels (a) and (c) are for soluble surfactant with
$k_{ads}C_{10}=12\,{\rm s}^{-1}$
and
$28\,{\rm s}^{-1}$
, respectively. Panels (b) and (d) are for an insoluble surfactant
$k_{ads}C_{10}=0\,{\rm s}^{-1}$
and differ only in the extent of collapse.
Appendix B.
In this appendix are examined in detail the gross characteristics of the main flow variables,
$h(x,t)$
,
$u_s(x,t)$
,
$Q(x,t)$
and
$\tau _w(x,t)$
, as the oscillation amplitude increases. The goal is to test the accuracy with which they may be roughly estimated from the linear solution.
First, the displacement of the interface from its equilibrium position is considered, as quantified by its extreme values
$h_{min}(x)-H(x)$
and
$h_{max}(x)-H(x)$
. These extrema are normalised with the local equilibrium value, forming the ratios
$\Delta h_{min}=[h_{min}(x)-H(x)]/H(x)$
and
$\Delta h_{max}=[h_{max}(x)-H(x)]/H(x)$
. The normalisation is motivated by the previous finding (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022) of a trivial solution (of the equations but not the boundary conditions) with purely radial velocities and displacement proportional to the equilibrium film thickness. Thus, the deviations of
$\Delta h_{min}$
and
$\Delta h_{max}$
from uniform values along the entire interface are indicative of the intensity of local shearing motion.
In the linear limit,
$\Delta h_{min}$
and
$\Delta h_{max}$
are proportional to the alveolar oscillation amplitude. Thus, the ratios
$\Delta h_{min}/a$
and
$\Delta h_{max}/a$
, plotted for
$a=0.01,\,0.05,\,0.10,\,0.15$
and
$0.20$
in figure 15(a), offer a first view of the deviations from linearity. It is readily observed that
$\Delta h_{min}/a$
and
$\Delta h_{max}/a$
vary strongly in the region from the rim and up to
$x\approx -0.5\;(\theta =60^\circ)$
, while they attain relatively uniform values deeper in the interior of the alveolus. This trend indicates that the shearing motion is more intense in the first
$30^\circ$
from the alveolar entrance.
(a) The maximum and minimum interface displacement normalised with the equilibrium film thickness and the oscillation amplitude. (b) The mean liquid film thickness normalised with the equilibrium thickness. Oscillation amplitudes
$a=0.01$
(dashed black),
$a=0.05$
(red),
$a=0.10$
(blue),
$a=0.15$
(green) and
$a=0.20$
(magenta).
In addition, it is observed from figure 15(a) that the symmetry around the equilibrium position (which holds for
$a=0.01$
) is gradually broken, with the liquid film thickening more during alveolar contraction and thinning less during expansion. This is further confirmed by figure 15(b), which displays the term
$\Delta \bar {h}(x)=[\bar {h}(x)-H(x)]/H(x)$
, with
$\bar {h}(x)$
the mean film thickness at location
$x$
. The increase in the mean film thickness is non-uniform close to the entrance, with the affected region increasing significantly with oscillation amplitude. The generally observed displacement of the liquid film thickness towards positive values is a purely geometrical effect because the variation in the interfacial surface area is not symmetric around its equilibrium value when the radius oscillates with finite amplitude.
Next, we examine the order of magnitude of the surface velocity,
$u_s(x,t)$
, volumetric flow rate,
$Q(x,t)$
and wall shear stress,
$\tau _w(x,t)$
. Figure 16(a,b,c) shows the maximum and minimum values attained by the above variables during an oscillation cycle, normalised with the oscillation amplitude. The amplitudes considered cover a wide range, starting from the linear limit,
$a=0.01$
, and going up to
$a=0.20$
, where extensive monolayer collapse occurs. It is observed that the results for the surface velocity are almost perfectly superimposed, indicating that the Marangoni stresses totally dominate the interfacial dynamics for all oscillation amplitudes. In contrast, the volumetric flow rate and the wall shear stress are also affected by capillary stresses and viscosity. (The effect of capillarity is clearly witnessed by the oscillatory variation of
$Q_{max},\,Q_{min}$
with
$x$
.) However, it is evident that the variations of
$Q_{max}/a,\,Q_{min}/a$
and
$\tau _{w,max}/a,\,\tau _{w,min}/a$
with amplitude are not dramatic, indicating that the role of oscillation amplitude is described by a regular expansion, i.e. higher-order terms provide small corrections. Thus, it appears that a convenient rule of thumb for rough estimates of magnitude is to use the linear results with any value of the amplitude
$a$
.
The maxima and minima of (a) surface velocity, (b) volumetric flow rate and (c) wall shear stress along the film for oscillation amplitudes
$a=0.01$
(dashed black),
$a=0.05$
(red),
$a=0.10$
(blue),
$a=0.15$
(green) and
$a=0.20$
(magenta).
Open access
