This section focuses on mathematical models, which are crucial for understanding the complex interaction between CLs and the ocular surface, analyzing tear film response to blinking and simulating CL fitting mechanics.

Different mathematical models have been subject to prior research and were dedicated to helping study the dynamics of tears and TBUT [12, 18]. Most tear film models are 1D single-layer models simplifying with consideration of the aqueous layer to be a Newtonian fluid [12, 18,19,20] and treating the tear film lipid layer as an insoluble surfactant monolayer [21, 22]. Researchers assumed the shape of the human cornea is negligible, and theoretical articles suggest using Cartesian coordinates on a flat substrate to develop models for the tear film and it is referred to as flat cornea approximation [19, 23]. Models were simplified to investigate important effects on tear dynamics such as evaporation and gravity over the open surface of the eye [19], osmosis across the corneal surface [12], Marangoni effects induced by varying lipid concentration [21], and complete and partial blinks [20, 24], among others.

Expanding on the dynamics of tears and TBUT from various mathematical models, research highlights the crucial parameter of tear film thickness, offering detailed analyses of the PCTF, PLTF, PoLTF, and the lipid layer [6, 7]. Contrasting earlier measurements, recent investigations estimate the human PCTF thickness at about 3 μm, although the thickness can vary [7, 25]. Post-blink, the tear film is influenced by surface tension gradients, and its TBUT is closely related to the thickness of the lipid layer, with factors like surface tension and evaporation playing roles [18, 26,27,28]. Meanwhile, Wong et al. [18] provided insights into the deposition process of the tear film, noting that the exposed eye section of the coating measures approximately 10 μm. Their model predicts film thickness and post-blink lipid spreading, showing that the film quickly thins at the edges and breaks when it becomes too thin, a process influenced by tear viscosity, initial thickness, and observed TBUT.

Specific studies focusing on the physical properties of the tear film highlighted key factors like viscosity [29] and surface tension [30], crucial for understanding tear film behavior and stability on the ocular surface. It found that healthy eyes typically exhibit a tear viscosity of about 6 mPa-s, in contrast to the higher average of 30 mPa-s in dry eye conditions, suggesting that tear film rheology could be significant in diagnosing and managing ocular issues, particularly for CL wearers [29]. Furthermore, the research revealed that tears from dry eye patients have increased surface tension compared to those from healthy eyes, contributing to a reduced TBUT [31]. This heightened surface tension, coupled with increased viscosity, leads to greater tear film instability and quicker tear film breakage, exacerbating dry eye symptoms [30]. These findings emphasize the complexity of tear film dynamics and the interaction of several factors in causing ocular discomfort and visual disturbances.

An in-depth analysis of tear fluid characteristics reveals its non-Newtonian nature due to the presence of molecules like proteins, lipids, electrolytes, and mucins [9, 32]. These components exhibit shear-thinning behavior [29], impacting tear film dynamics when lipids are removed [33]. Some studies have neglected the influence of the corneal curvature, assuming a spherical substrate shape [34], while others explored cylindrical or prolate spheroid geometries [35, 36]. While the prolate spheroid approximates the human cornea, research suggests that corneal shape has minimal impact on tear film thinning rates, often leading to the assumption of a flat cornea in computational models of tear film dynamics [23]. Understanding tear fluid properties and substrate curvature helps refine tear film models.

The exploration of tear film dynamics begins with an examination of parameters related to the ocular surface, emphasizing mathematical models that detail tear film formation and relaxation during blinking, as depicted in Fig. 2 [19]. These models, focusing on tear film thickness [7], delve into the evolution of the aqueous layer and consider factors like evaporation and heat transfer [19]. Typically, simulations begin with an initial condition, assuming uniform tear film deposition except for menisci near the eyelids. While many models simplify the lipid layer by assuming a stress-free upper surface [20, 22, 24], evidence suggests its role in particle movement [21, 37, 38]. Theoretical studies explore mathematical models with different blinking characteristics, incorporating upper lid movement during the opening phase of the eye [18, 22]. These models introduce fluxes to estimate tear supply, concluding that no-flux conditions fail to provide adequate coverage, necessitating tear fluid flux from the eyelid. Additionally, the impact of the lipid layer highlights altering the distribution of the film, influencing film height between blinks [18, 37]. These models also introduce critical concepts like the stress-free limit (SFL) and uniform stretching limit (USL), marking milestones in tear film research.

Schematic diagram from a mathematical viewpoint: (a) PCTF; (b) PLTF. (Parameters: X(t) – position as a function of time t; L – half-width palpebral fissure; h0 – initial tear meniscus height from both eyelids; Hcl – thickness of CL; D – PoLTF thickness; hPCTF, hPLTF – PCTF thickness and PLTF thickness respectively as a function of time)

Another tear film model was devised to investigate tear film evolution across multiple blink cycles, focusing on the opening and closing of the eye. Lid movements were characterized using two approaches: sinusoidal motion [20] and realistic blinking [24]. Braun and King-Smith employed a sinusoidal blink model [20], defining each cycle as a period and identifying periodic solutions for both complete and partial blinks. This model revealed distinctions between the superior and inferior tear film layers and replicated in vivo tear film observations during partial blinks, offering insights into defining complete blinks based on fluid dynamics principles [20].

In contrast, the realistic blink model incorporated the entire blink cycle, including lid opening, open-eye duration, and closure, based on actual lid motion data [13, 39]. This approach introduced flux boundary conditions accounting for lacrimal gland supply and punctal drainage. Results indicated thicker tear film near the moving end during lid opening and closing, with thinning near the ends during stationary, fully open phases of the eye, providing a more comprehensive view of tear film dynamics [24]. Recent advances in imaging and mathematical models of tear film dynamics during blinking have deepened our understanding, with increased use of simultaneous imaging and improved OCT instruments promising further insights [40].

Mathematical models enhancing our understanding of tear film dynamics have concentrated on the lipid layer, particularly modeled as a polar lipid monolayer [22, 37]. These models demonstrate how concentration gradients induce Marangoni flow and affect evaporation rates within the tear film lipid layer, which is sensitive to variations in pressure, temperature, and surfactant concentration [41]. Although effective in capturing TBUT due to evaporation, these models struggle with identifying increased evaporation rates influenced by surfactant concentration. In scenarios considering both fully open and half-closed eye states, they compute high lipid concentrations near the lower lid during interblink, propelling the lipid upward and consequently dragging the aqueous tear fluid. When the eye half-blinks, this concentration peaks at the center of the eye, causing rapid tear film thinning. Another model examines the evolution of tear film thickness and lipid concentration during blinking [21], revealing that higher lipid concentrations amplify the Marangoni effect, driving lipids toward the upper lid. This model, framed as a coupled partial differential equation, shows that the presence of lipids not only thickens the tear film due to increased fluid flow but also results in a non-uniform lipid distribution across the tear film.

Tear film dynamics research has undergone significant evolution, introducing a model with a lipid reservoir continuously supplying lipids to the system and altering boundary conditions to control the impact of the lipid flow on tear film evolution [42]. Simultaneously, studies have explored rapid tear film thinning linked to uneven lipid layer distribution, emphasizing the correlation between a healthier, more uniform lipid layer and an extended TBUT due to lower surface tension [43]. Furthermore, investigations have utilized fluorescein imaging to simulate tear film thinning and solute transport, aligning simulated fluorescein intensity with in-vivo observations to differentiate between evaporative and tangential flow-driven tear thinning mechanisms. Prior research also qualitatively observed tear film thickness distributions and a drop in polar lipid content near the lids during blinking [21]. Additional model variations incorporated a dilute surfactant model and a thick aqueous layer with large menisci, revealing evidence of substantial lipid remnants after the upstroke of the blink cycle and the potential for a significant boundary thickness to facilitate tear film development during a full upstroke [38].

Briefly, the dynamics of the tear film are crucial for ocular health, and mathematical models provide invaluable insights into the layered complexity of the tear film. In the next section, we will explore computational solutions to the tear film dynamics in the presence of CLs, further advancing the knowledge in this field.

The interaction between the tear film and CLs is an intricate aspect of ocular physiology, as depicted in Fig. 2b, with implications for lens comfort and visual acuity. Investigating the dynamics of tear film behavior in the presence of CLs provides insights into the optimal design. Central to this exploration is the mathematical modeling of tear film behavior. Expanding the tear film models, recent advancements have incorporated specific parameters pertaining to CL wearers. Utilizing a lubrication theory-based approach, these models adeptly describe the dynamics of tear film in the context of blinking and CL wear [44,45,46,47,48]. It is worth noting that in the mentioned studies, the tear film is bifurcated into two distinct layers: the PLTF, which is the fluid layer sandwiched between the CL and the external environment, and the PoLTF, situated between the CL and the corneal surface.

Hayashi and Fatt [44] used lubrication theory to investigate tear exchange caused by the compression of a soft CL by the eyelid against the cornea, finding that each blink leads to an estimated 10–20% tear exchange with typical film thicknesses of 8 to 10 microns, highlighting the importance of blinking for maintaining tear film balance in CL wearers and the need for understanding tear dynamics to improve CL design and guidelines. Building on these findings, subsequent research delved into tear film dynamics with CLs, considering factors like lens thickness, permeability, gravitational effects, and slip models at the fluid-lens interface [45]. A complex mathematical model, derived by applying a lubrication approximation to hydrodynamic motion equations and considering the porous layer of the tear film, was developed to study the post-blink film evolution, revealing that increased lens thickness, permeability, and slip could accelerate film thinning, although these changes have minimal effect under standard CL conditions.

Dunn et al. [46] investigated the impact of blinking on tear film dynamics with soft hydrogel CLs, observing that blinking can either partition the tear film or fully integrate it into the CL, leading to relative sliding between the lens, corneal epithelium, and eyelid wiper. Their numerical model shed light on the pressures and sliding speeds involved, emphasizing that eyelid-lens interaction predominantly occurs in a hydrodynamic regime and is critical for understanding the lubrication behavior of CLs, particularly regarding ocular sliding, loading, and the potential for surface damage due to shear stress [2]. Building on these findings, Talbott et al. [47] explored the effects of evaporation on the PLTF with permeable CLs, noting how evaporation reduces PLTF thickness and leads to fluid loss through the lens. They employed lubrication theory to formulate an equation representing PLTF thickness, accounting for evaporation, thermal transfer, and capillary action. Their study compared comprehensive and simplified models, offering insights into fluid loss due to evaporation and contributing to a deeper understanding of the fluid dynamics involved in CL wear.

Anderson et al. [48] furthered the understanding of tear film dynamics with CLs by examining the partitioning of the PCTF similar to prior research, with thicknesses ranging from 1 to 5 µm, in contrast to the considerably thicker CLs (50–400 µm). They noted that CLs are subject to forces in both horizontal and vertical directions during blinking, with recent studies focusing more on vertical movement. Chauhan and Radke [49] evaluated this vertical motion using an innovative method based on mechanical force balance, considering forces from the eyelids, gravity, elasticity, and viscosity, and integrating parameters like lens attributes and tear film thickness variations. They discovered that the downward movement of the lens during a blink is 2–3 times greater than during the interblink phase, indicating that current testing methods may overlook significant aspects of lens movement. This research emphasizes the need for more comprehensive experimental approaches in understanding CL behavior and tear film dynamics.

Maki and Ross [50] introduced a novel method to calculate the suction pressure under a soft CL, focusing on how the lens deforms under the combined forces of the tear film and eyelid blink. Their findings revealed that with a consistent eye shape, the center of the lens experiences more suction pressure as the curvature radius of the lens increases, while peripheral pressure decreases, and negative pressure in the transition zone increases for larger radii. Building on this, one research [51] examined the impact of CL design on ocular health and how blinking affects lens adaptation, particularly how the lens attempts to regain its shape and generates suction in the PoLTF, influencing tear fluid movement and potential fluid exchange at the lens edge. Another study [52] employed a variational method to assess elastic stresses in CLs and their associated suction pressure, providing solutions to the Euler–Lagrange equation for lenses with consistent thickness, although challenges arise with variable lens thickness. These studies advance the understanding of the forces at play in tear film dynamics and CL behavior.

In summary, tear film dynamics with CLs merge biomechanics and fluid dynamics, with foundational research offering insights into tear film behavior, lens design impacts on ocular health, and underlying mathematical models. The next section will focus on the methods and challenges in solving these models, enhancing our understanding of the topic.

Tear film dynamics involve creating non-linear partial differential equations with appropriate conditions, and MATLAB is a commonly used platform for solving these models [53]. MATLAB simplifies complex data analysis, offering a programming language for numerical computations and mathematical tasks with functions for matrices, algorithms, and user interfaces. Researchers employ MATLAB to solve mathematical models of tear film dynamics, often using the finite difference method [54], which is straightforward but may require a high number of grid points, making it time-consuming [19, 20, 22, 37, 38]. Despite its simplicity, this approach demands additional studies on stability and accuracy.

Alternatively, the Chebyshev spectral collocation method provides a more advanced solution [55]. It utilizes non-symmetric mapping to minimize point spacing, transforming equations into a time-dependent system of differential algebraic equations. This method further enhances accuracy and computation speed through a modified non-symmetric mapping and two input parameters, reducing errors, particularly in higher-order derivatives [56]. These techniques offer researchers powerful tools to analyze tear film dynamics effectively, balancing ease of use and computational efficiency [24, 48].

Concisely, the analysis of tear film dynamics using computational methods emphasizes the importance of mathematical models in predicting tear film behavior, leading to a subsequent section on the mechanical interplay between the cornea, CL, and ocular surface.

The simulation of the computational mechanical models gives an understanding of the impact of the fitting of CL shape over the eye and their corneal pressure and friction. There are not many available models of this type in a Multiphysics approach due to the disparity in in-vivo measurement parameters and real characteristics of complex fluid dynamics of tear film [2, 46, 57,