Pattern formation, structure and functionalities of wrinkled liquid crystal surfaces: A soft matter biomimicry platform

This review presents an integrated theoretical and computational characterization and analysis of surface pattern formation in chiral and achiral liquid crystal self-assembly and the mechanical/optical/tribological/tissue engineering surface functionalities that emerge from various wrinkling processes. Strategies to target surface patterns include linear, non-linear, multidirectional and multiscale wrinkling phenomena. The focus of the review is to show the unique surface structure-functionalities that emerge from anisotropic liquid crystal soft matter, eliminating or reducing the need of aggressive solvents, extreme pressure/temperature conditions, erosion and other surface morphing approaches. The surface pattern formation theoretical-modelling- computational results are then connected and validated with actual biological surfaces that are considered solid liquid crystal analogues, such as exocuticles of insects, fish scales, and flowers. A unique feature of the in silico surface pattern formation platform used throughout this review is the generalized liquid crystal shape equation that includes surface anchoring elasticity, membrane elasticity, and stress loads from liquid crystals orientation gradients. Clear characterization of surface shapes, curvatures, roughness, that are behind surface functionalities are introduced and applied to strengthen validation of predictions with actual nature’s surfaces. Wrinkling scaling laws, and the dependence of material properties on morphing mechanisms are elucidated. The predictions capture very well the two-scale wrinkling patterns in tulips, wrinkling gradients that display water sensor capabilities, egg carton shapes in rose petals and their potential for cell alignment, and the ability to create surface roughness with targeted kurtosis and skewness to control and optimize friction and tribological functionalities. The results are summarized in terms of surface geometry (open or closed) mechanisms and phenomena (anchoring, membrane elasticity), material properties (anchoring coefficients, membrane bending modulus, Frank elasticity), wrinkling scales and scaling laws (amplitude, wave-lengths, skewness, kurtosis) and functionalities (optical iridescence, friction, wettability, structural color, curvature-driven cell alignment and differentiation). Taken together, the range of surface geometries and surface functionalities captured by the liquid crystal biomimetic in silico platform provides a foundation for future experimental green manufacturing pathways based on anisotropic soft matter.


Introduction
Liquid crystallinity is a state of matter between solid and liquid when rod-like molecules align along a preferred macroscopic orientation (characterized by a unit vector or director n), exhibiting non-Newtonian viscous fluidity and anisotropy (Rey and Tsuji, 1998;Rey and Denn, 2002;Rey, 2009). Solid liquid crystal analogues are frozen liquid crystals due to external stimuli such as the presence of crosslinkers (Liang et al., 2021), and laser or plasmatriggered polymerization Kameche et al., 2020). Solid liquid crystal analogues are ubiquitous in biological systems (Mazur et al., 1982). For example, chitin fibres with a helical structure form the exoskeleton of shimmering beetles (Sharma et al., 2009;Mitov, 2017), responsible for the structural colour. A similar structure is also found in mantis shrimps Weaver et al. (2012); Patek et al. (2004) and Homarus americanus (Raabe et al., 2007;Nikolov et al., 2015), which plays a significant role in their exceptional mechanical endurance. Biological plywood and human compact bones are other typical liquid crystal analogues, whose cellulose and collagen fibres align in a chiral liquid crystalline-like hierarchy, known as the Bouligand architecture (Cowin, 2004;Lagerwall et al., 2014;Gutierrez and Rey, 2017). Cellulose-based liquid crystals have various applications (Habibi et al., 2010;Li et al., 2021) such as serving as a steam engine driven by humidity , biofibres' surface morphology sensor (Aguirre et al., 2016), piezoelectric sensor (Csoka et al., 2012;Rajala et al., 2016), water-based lubricant (Noroozi et al., 2014;Shariatzadeh and Grecov, 2019;Zakani et al., 2022), inducing long-range order in a semiconducting polymer (Risteen et al., 2017). The unique spatiotemporal and structural colours Saraiva et al., 2020;Silva et al., 2021) makes celluloses a promising optical material (Lagerwall et al., 2014;Ličen et al., 2016;Anyfantakis et al., 2020). Other examples such as DNA both in vivo (Bouligand, 1972) and in vitro (Chow et al., 2010), viruses Ilca et al. (2019), spider silk (Saravanan, 2006;Kluge et al., 2008), fibroblasts and osteoblasts (Kemkemer et al., 2000a;Kemkemer et al., 2000b) are also natural liquid crystal analogues in biological system. Theory, computation Murugesan and Rey, 2010b;Murugesan and Rey, 2010c;Abou-Dakka et al., 2012;Khadem and Rey, 2019a;Khadem and Rey, 2019b;Khadem and Rey, 2021), and modelling (De Luca and Rey, 2003;Murugesan and Rey, 2010a;Gutierrez and Rey, 2014a;Gutierrez and Rey, 2014b;Murugesan et al., 2015;Aguilar Gutierrez and Rey, 2016;Gutierrez and Rey, 2016) of these systems are evolving and active fields seeking to link biological processes and liquid crystal self-assembly and liquid crystal self-organization (Park et al., 2021) to eventually contribute to biomimetic engineering of functional surfaces (Monaenkova et al., 2012;Lehnert et al., 2013;Tsai et al., 2014). The approach to link natures' shapes and bionsipired functional surfaces is liquid crystal capillary science: Various classes of surface corrugations were observed in liquid crystal analogues (Tan et al., 2020) including creasing, folding, wrinkling (Yoo and Lee, 2005;Wu et al., 2007;Cutolo et al., 2020), and ridges (Jin, 2014;Wang and Zhao, 2015). Surface creasing is mostly studied in swelling-induced hydrogels (Trujillo et al., 2008;Zhou et al., 2017;Dortdivanlioglu et al., 2021), elastomers (Hong et al., 2009;Chen et al., 2012;Liu et al., 2019;Guan et al., 2022) and growth-driven soft tissues (Jin et al., 2011). Folding emerges in biological membrane system, such as in the Golgi apparatus (Norlén, 2001), mitochondrion membranes (Perkins and Frey, 2000;Joubert and Puff, 2021) and cortical folding (Mangin et al., 2004;Striedter et al., 2015;Garcia et al., 2018). The morphologies of biological membranes determine their functionalities. For example, the smooth endoplasmic reticulum sythesizes lipids and phospholipids (Cartwright and Higgins, 2001;Wilkinson et al., 2002), while folded rough endoplasmic reticulum is responsible for the secreted proteins (Lodish et al., 1983;Mandon et al., 2013). The surface corrugations usually perform a multiscale structure. Compared to creasing and folding, surface wrinkling has a smoother structure. Creasing, folding and wrinkling patterns can interconvert. For example, thermo-responsive wrinkling patterns can transform into foldings and can be applied to soft actuators (Tokudome et al., 2016). For load-driven surface corrugations, the creasing pattern appears at a lower critical compressive loading compared with the wrinkling pattern, which have both been validated experimentally (Gent and Cho, 1999;Cai et al., 2012) and theoretically (Karpitschka et al., 2017;Ciarletta, 2018;Ciarletta and Truskinovsky, 2019). Surface pattern formation, characterization and functionality are formulated using differential geometry (Koenderink and Van Doorn, 1992;Kornev and Shugai, 1998;Kornev et al., 1999;Gutierrez and Rey, 2018;Zhang et al., 2018;Wang et al., 2020b;Alimov et al., 2021). A local surface patch has two principal surface curvatures κ 1 and κ 2 (κ 1 < κ 2 ), which can be used to describe surface geometry, as in the mean curvature H = (κ 1 + κ 2 )/2 and the Gaussian curvature K = κ 1 κ 2 . The signs of (H, K) are used to classify surfaces such as cylinders, domes, and saddles. These geometric quantities (H, K) comingle shape and length scales since H and K have dimensions of the reciprocal length but another measure can be used that avoids this commingling of shape and length scale. Hence it is useful to use the dimensional Casorati curvedness C) and the normalized, dimensionless shape parameter S), developed in the previous literature (Koenderink and Van Doorn, 1992;Gutierrez and Rey, 2018;Wang et al., 2020b). The curvedness characterizes the amplitude of wrinkling, while the shape parameter is a dimensionless scalar that differentiates spherical, cylindrical and saddles patches. In this manuscript, we focus on 1D and 2D open surface corrugations. A 1D surface corrugation (for example, Figure 1C) shows a wrinkling profile along a single direction (x or y) and has a vanishing Gaussian curvature. A 2D surface corrugation ( Figure 1F) depends on two orthogonal directions (x and y). Complex synthetic and biological surfaces have a spatial distribution of amplitudes and shapes and hence distribution functions are required. A global statistical method can be applied to study surface roughness (Tayebi and Polycarpou, 2004;Hansson and Hansson, 2011;Sedlaček et al., 2012). The local curvatures and global geometric statistics determine the surface functionalities. For example, a cytoskeleton organization can be topology-induced on a local concave (S = −1), convex (S = +1), and hyperboloid surface (S = 0) (Yang et al., 2022). Torn plastic sheets and leave's photosynthesis are also geometry-driven wrinkling process (Nath et al., 2003;Sharon et al., 2004;Sharon et al., 2007). Other examples of curvaturedependent applications in tissue engineering (Rumpler et al., 2008;Blanquer et al., 2017) can be found in regulating stem cell migration and differentiation (Werner et al., 2017), the behaviours of cellular growth  and early osteogenesis Ozdemir et al. (2013), bone tissue regeneration (Zadpoor, 2015) and scaffolds with various geometry (Bidan et al., 2013).
There are multiple ways to generate surface corrugation, such as the presence of plasma (Nania et al., 2015;Cheng et al., 2016;Cheng et al., 2017) or laser (Hu et al., 2017;Wang J. et al., 2019), electric (Lin et al., 2020) or magnetic field Wu and Destrade (2021); Tokarev et al. (2014), and stress field (Rodriguez-Hernandez, 2015). The change of humidity , pH González-Henríquez et al. (2018), temperature , polymerization (Bergmann et al., 2020;Junisu et al., 2022) and drying process Izawa (2017) also disturb surface wrinkling. Surface corrugations are widely used in bioinspired multifunctional materials, such as the flexoelectric materials (Castles et al., 2011;Herrera-Valencia and Rey, 2014;Rey et al., 2014b;Herrera Valencia and Rey, 2018;Zhang et al., 2021) and curvature-induced magnetochiral materials (Hertel, 2013;Sloika et al., 2014). A liquid crystal surface coated with surfactants shows elasticity, whose classical theory is first proposed by Helfrich (Helfrich, 1973), and further completed by Ou-Yang (Zhong-Can and Helfrich, 1987). The membrane elastic model successfully predicted the biconcave discs observed in human red blood cells and was validated in the literature (Seifert et al., 1991;Kuzman et al., 2004;Mu et al., 2020;Martínez-Balbuena et al., 2021). Theoretical research of surface wrinkling patterns by Helfrich elastic theory is studied in (Baumgart et al., 2005;Bian et al., 2020). However, the elastic model of a system with open free edges is problematic. Tu (Tu, 2010;Zhan-Chun, 2013) proved that if line tension is finite, an open membrane cannot have constant non-zero mean curvature (Tu and Ou-Yang, 2014). More numerical solutions to surfaces' shapes by using a pure elastic model are studied in (Umeda et al., 2005;Li et al., 2013). However, all the models discussed above require an external force field, such as an electric (Lin et al., 2020) or magnetic field Wu and Destrade (2021), or the presence of mechanical stress (Rodriguez-Hernandez, 2015). The anchoring model (Rey, 2010;Rey and Herrera-Valencia, 2012;Rey et al., 2014a) on the other hand, enables the discussion of selfassembled structure, which is a unique architecture observed in liquid crystal. The anchoring model (Rey, 2007;Rey and Murugesan, 2011;Rey et al., 2013) is a prevailing method when multiscale wrinkling patterns are formed. The anchoring model is based on the anisotropic nature of the interfacial tension of liquid crystals (Rey, 1999a;Rey, 2000a;Rey, 2000c;Cheong et al., 2001;Rey, 2001b;Cheong and Rey, 2002a;Cheong and Rey, 2002b;Cheong and Rey, 2002c;Zhou et al., 2002;Rey, 2007;Kim et al., 2016). The total interfacial (surface) tension involving a liquid crystal of orientation n (n = −n, n 2 = 1) at a surface with outward unit normal k is: γ = γ 0 + μ 2 (n ·k) 2 + μ 4 (n ·k) 4 + . . .. The presence of liquid crystalline orientation n) and anisotropic surface tension μ 2 (n ·k) 2 + μ 4 (n ·k) 4 + . . . creates a new pathway to surface pattern formation models and predictions based on liquid crystalline self-assembly and self-organization. For example, capillary pressures can include dilation, tilting and orientation contributions, enriching the number and types of forces and torques that can be used to create wrinkling and folding patterns, close to those observed in nature, likewise tangential forces due to anchoring and orientation gradients generate new types of Marangoni flow beside the usual thermo-and chemo-capillarity (Rey, 1999c). Generalizations to liquid crystal-membrane models, that combine anchoring and membrane elasticity yield shape diagrams that naturally incorporate the three shaping modes of bending, creasing, and folding (Yoo and Lee, 2005;Wu et al., 2007;Jin, 2014;Wang and Zhao, 2015;Cutolo et al., 2020) and point to a rich unexplored area of surface pattern formation in soft matter. In this review we summarize key principles and most significant applications of liquid crystal capillarity applied to synthetic and biological surface pattern formation, building on our previous contributions (Rey, 1999a;Rey, 2000a;Rey, 2000c;Cheong et al., 2001;Rey, 2001b;Cheong and Rey, 2002a;Cheong and Rey, 2002b;Cheong and Rey, 2002c;Rey, 2007) and shedding light on the origin of scaling laws, material properties, and optical/sensor/actuator/wetting/ tribological functionalities. Some of the issues of key phenomena, processes, mechanisms, and functionalities presented in this review are.
1. What are the new capillary pressure contributions to shape generation arising from liquid crystal anisotropic interfacial tension? 2. What are the scaling laws relating orientation n) length scales and anchoring coefficients W), and wrinkling amplitudes and wrinkling wavelengths? 3. What are the new phenomena arising from higher order (quartic as opposed to just quadratic) formulations of the anchoring energy, including multiple wrinkling scales and amplitudes? 4. What is the role and impact of surface orientation gradients (∇ surface n) on shape generation? Is there a connection between basic orientation deformation modes (known as splay, bend) and mean surface curvature H)? 5. What is needed to generate biaxial wrinkling, such as those found in egg cartons and their variants? 6. What are typical surface geometry statistics (important in tribology, friction, and wetting), including kurtosis and skewness, emerge from different wrinkling patterns? 7. What are distinct surface functionalities, such as optical, tribology, wetting, that are captured by the liquid crystal anchoring-driven wrinkling? 8. What are the new surface pattern formation modes that emerge by combining liquid crystal anchoring and membrane bending elasticity?
Answering these questions provides the framework to develop bioinspired functional surfaces via the liquid crystal biomimetic platform. Figure 1 presents a flowchart of the organization and flow of information in this review. The anchoring-driven surface wrinkling patterns appear in both open surfaces (such as tulip petals) and closed surfaces (such as drops or tactoids). In Section 1 we introduce the background and motivation for studying the surface of liquid crystal analogues found in nature. In Section 2.1 we review the 1D open surface corrugation governed by the pure anchoring model and the application to friction. In Section 2.2, we review the influence of the helix pitch of cholesteric liquid crystal analogues on the geometry and the application in diffraction gratings. The membrane elastic model is introduced in Section 2.3 for obtaining a more complex wrinkling pattern, with application in surface wettability. In Section 3, we focus on the 2D open surface corrugations. In Section 4, we review the wrinkling profile of the liquid crystal droplet with an enclosing elastic membrane. The closed droplet morphology reviewed in Section 4 shows a completely different wrinkling mechanism compared to the open surface discussed in Section 2 and Section 3. The bulk-surface interaction plays a significant role in the wrinkling formation of a closed droplet, which can be neglected for an open surface. In this review, wrinkling is emphasized but creasing and folding are also discussed.
2 1 D open surface wrinkling: Cylindrical patterns 2.1 Pure anchoring model with application to friction The multiscale surface wrinkling is a unique structure observed in many plants (Lin et al., 2017;Schauer et al., 2017), for example, the Queen of the Night tulip (Oh et al., 2019), Yunnan rhododendron, daisy Ursinia calendulifolia (Antoniou Kourounioti et al., 2013) and Hibiscus trionum (Vignolini et al., 2015;Chen et al., 2021). The tulip Queen of the Night exhibits glossy and iridescent colour due to the scattering and diffraction effect ( Figure 2B). The special optical effect for those petals plays an important role to attract bees Whitney et al. (2009), further benefiting the pollination van der Kooi et al. (2015). The important point here is that the surface is composed of larger and smaller features that can be approximated by two-scale wrinkling. Figure 2 presents the biological example and numerical simulation for 1D open surface wrinkling. Figure 2A shows the AFM image of the cylindrical surface corrugations observed on Tulip Queen of the Night. Figure 2B is the schematic of how the wrinkling geometry affects the glossy and iridescent colour of tulip petals by light scattering and diffraction.
The anchoring model is developed to describe the multiscale cylindrical surface wrinkling profile and the pattern formation mechanism (Rofouie et al., 2014;Rofouie et al., 2015a;Rofouie et al., 2015b;Rofouie et al., 2018;Wang Z. et al., 2019;Wang et al., 2020a). A cylindrical surface means 1D wrinkling as in a simple corrugation, with one-directional wave-vector. As mentioned above, the anisotropic surface energy density γ is described by the Rapini-Papoular equation, which regards γ as a function of the director n of liquid crystal, and the surface unit normal k Stewart (2019) (2007): where γ 0 is the isotropic surface tension, μ 2 * is the dimensionless quadratic anchoring coefficient, and μ 4 * is the dimensionless quartic anchoring coefficient. Eq. 1 implies that by reversing n → −n or k → −k, the surface energy density is unaffected. The contribution from higher order terms can be neglected due to −1 ≤ n ·k = cos θ ≤ 1, where θ is the angle between surface normal k and director n, as shown in Figure 3. Planar anchoring ( Figure 3A) and homeotropic anchoring ( Figure 3B) correspond to θ = π/2 and θ = 0, respectively. Figure 3D shows the normalized surface tension reduction as a function of θ by assigning different μ 4 /μ 2 values. The maximum surface tension reduction occurs at a different position, implying a complex wrinkling phenomenon dependent on μ 4 /μ 2 . The stable surface wrinkling profile is controlled by the momentum balance equation along k-direction (Rey, 1997a;Rey, 1997b;Rey, 1999a;Rey, 1999c;Rey, 2001a;Rey, 2001b;Wang Z. et al., 2019): where p c is the pressure difference between the liquid crystal phase and the isotropic phase, f b is the bulk elastic energy density, ξ is the Cahn-Hoffman capillary vector (Hoffman and Cahn, 1972;Cahn and Hoffman, 1998). The capillary vector ξ is a useful tool to model surface pattern formation and can be seen as the sum of a Lagrange multiplier or constraint (normal component) and the Lagrangian derivative (tangential component): where F s is the total surface free energy obtained from the area integral of the interfacial tension. An open surface has a vanishing stress jump (Wang et al., 2020a), which results in the governing equation (Wang et al., 2022) 0 Eq. 4 captures the wrinkling patterns on the tulip petals (Rofouie et al., 2014;2015b;a, 2018), which are the cancellation effects between the dilation pressure, director pressure and rotation pressure. The director pressure depends on the director field n, which is the driving force of the wrinkling profile. To eliminate the director pressure, a surface patch induces a local rotation (rotation pressure) and local dilation (dilation pressure), which changes the geometry by affecting k, and is performed as self-assembly. Figure 2C is the anchoring parameter space μ 4 − μ 2 , superposed with the computed surface wrinkling profiles h(x). The wrinkling parametric space is characterized by a resonant line (thick red line) and a damping line (thin black line). Maximal two-scale wrinkling occurs above and below the resonant line. Two-scale wrinkling only arises in quartic anchoring models when the anchoring coefficients have opposite signs (−+) or (+−). Figure Figure2d) shows representative two-scale wrinkling profiles and defines the large (h 0 )and small (h 2 ) scales of the amplitudes. The solution of the shape equation and analysis (Rofouie et al., 2018) reveals the presence of novel quadratic scaling laws for two-scale h 0 , h 2 wrinkling: h 0 h 2 μ 2 2μ 4 2 1 + μ 2 2μ 4 2 (5) whose critical value h 0 /h 2 = 2 is easily accessible by the anchoring coefficients. Four fundamental wrinkling modes: two waves, flat hill/ valley, four waves and periodic double are summarized in Figure 1C.
In partial summary these results clearly show that multiscale wrinkling emerging from higher order anchoring energies is a simple manifestation of the coupling between liquid crystal orientation n) and local surface geometry k) and the scaling law and up-down symmetry is hence fertile in surface functionalities.
Having described and analyzed 1D two-scale wrinkling due to higher order anchoring, next we discuss functionalities arising from rough, non-planar surfaces that we are able to produce in silico using the shape Equation 4. The surface wrinkling is often studied along with the roughness characterization due to its effects on surface friction, restitution, and adhesion Persson et al. (2004); Emmens (1988); Li et al. (2020a). Two significant roughness parameters are the skewness (Ssk) and kurtosis (Sku), corresponding to the third-order and fourthorder moment information of wrinkling profile distribution function (Tayebi and Polycarpou, 2004). Skewness and kurtosis are dimensionless parameters that are scaled by the third and fourth order of standard deviation. A normal distribution results in a vanishing skewness and Sku = 3. The higher the absolute value of skewness, the larger the deviation from the normal distribution. Positive skewness implies a surface with tall peaks, and negative skewness implies a surface with low valleys (see Figure 4B). The kurtosis portrays how much sharp peaks or valleys contribute to the surface profile (see Figure 4C). The skewness-kurtosis plot ( Figure 4D) can be used to study the roughness of biological surface topography (Wainwright et al., 2017), thermal-hydraulic performance Liao et al. Frontiers in Soft Matter frontiersin.org 06 (2020,2021), normal contact stiffness (Belhadjamor et al., 2020) and other processes Ba et al., 2021). Figure 4a) shows representative results of skewness and kurtosis using Eq. 4. The wrinkling pattern shows rich skewness and kurtosis information when the two anchoring coefficients have opposite sign. A bioinspired application is found on fish scalps; Figure 4E shows the cylindrical wrinkling profile of spotted tinselfish, whose skewness and kurtosis data are evaluated in (Wainwright et al., 2017). Most of the surfaces are approximated close to a normal distribution, showing an almost vanishing skewness, and such surface can be easily obtained by setting μ 2 = μ 4 in pure anchoring model (see Figure 4A). Figure 4D shows that the anchoring model produces a wide range of skewness and kurtosis, with maximum values max (Rku) ≈ 3 and max (Rsk) ≈ 1. Positive skewness implies a surface with sharp peaks, which is mostly obtained by setting μ 4 > 0 and μ 2 < 0. Such surface is observed in triggerfish (Rsk = 0.42), armored catfish (Rsk = 0.13), and trout with mucus (Rsk = 0.15). On the contrary, negative skewness (μ 4 < 0 and μ 2 > 0) is observed in squirrelfish (Rsk = −0.11) and hammerhead shark (Rsk = −0.14). Most biological surfaces also show a kurtosis smaller but close to 3, implying a surface lacking of extreme peaks or valleys, for example, bonefish (Rku = 2.9), squirrelfish (Rku = 3), armored catfish (Rku = 2.8). Some extreme values of kurtosis can also be found in longnose butterflies (Rku = 4.1) and red maple leaves (Rku = 4.3) (Wainwright et al., 2017).
In partial summary, moving in the anchoring parameter space generates ( Figure 2C) a kurtosis-skewness bilobal loop (Figure 4d)) that captures a wide range of functional biological surfaces, including butterflies, maple leaves, and tinsel fish. The experiment validation can be found in (Maeda, 1999;Kirkwood and Fuller, 2009;Fernandes et al., 2013).

Surface with varying amplitude and optical applications
The most abundant examples of solid liquid crystal analogues in nature display the Bouligand architecture, or cholesteric liquid crystal (Patek et al., 2004;Sharma et al., 2009;Weaver et al., 2012;Sharma V. et al., 2014;Mitov, 2017). The helix pitch of cholesteric liquid crystal P 0 is the distance under which the chiral director field n undergoes a full rotation: n(x) = n (x + P 0 ). In Section 2.1, P 0 was considered as a constant, which in nature, could be sometimes unrealistic. The humidity of the environment changes the helix pitch, which further affects the surface profile and its optical properties. For example, beetles Dynastes hercules, Coptocycla, Tmesisternus isabellae (Hinton and Jarman, 1972;Hadley, 1979;Liu et al., 2009) and vascular plant Selaginella willdenowii (Lee and Lowry, 1975) can adjust their structural colour to the humidity of environment.
The Frank energy density is dependent to the Frank constant K (Priest, 1973;Straley, 1973), the helix pitch P 0 and the director field n: (Rofouie et al., 2015a). Including f b and solving Eq. 2 yields the numerical solutions to the wrinkling profile with varying helix pitch P 0 x). Figure 5 shows representative numerical simulations of surface wrinkling profiles with varying amplitude. Figure 5a) demonstrates how the following four significant factors appearing in our shape equation manipulate the wrinkling profile (Rofouie et al., 2015a): 1) the dimensionless quadratic anchoring coefficient W/γ 0 (equivalent to 2μ* in Equation 1, 2) the elasto-capillary length scale K/γ 0 , 3) the micron scale pitch in dry state P dry , and 4) the hydration-driven pitch gradient (P wet − P dry )/L, respectively. Figure 5A to Figure 5D verifies that by increasing all the four factors, the amplitude is enhanced. The surface amplitude doubles if the anchoring coefficient is also doubled. However, multiplying K/γ 0 by  Figure 5C shows that higher P dry results in a larger periodicity. And Figure 5D shows the both the amplitude and periodicity increase. They are dependent on the x-axis due to the gradient effect.
The periodicity shown in Figure 5A to d is around 1,000 nm, which is the same order as the wavelength of visible light. If the visible light is the incident light, the iridescence and colours appear due to diffraction grating ) (a schematic is previously shown in Figure 2B). Figure 5F presents the scattering patterns for a surface with varying pitch ( Figure 5E) with blue incident light (wavelength: 480 nm). The scattering and diffraction effect of surface wrinkling on beetles' exoskeletons cause iridescent colours (Sharma et al., 2009). The sensitivity of structural colour to humidity of beetles (Kim et al., 2010;Sun et al., 2021) is directly generated by varying pitch. The simulation result is validated by (Terris et al., 1992), as shown in (Rofouie et al., 2015a).
In partial summary, the general liquid crystal shape Equation 4 contains the mechanisms to generate variable spatial wrinkling due to spatial gradients in water content. The wrinkling can be used as water sensors when coupled to optical responses. This is another example of liquid crystal surface functionality that bio-mimics structural colour in nature.

Elastic anchoring model with application in hydrophobicity
The Rapini-Papoular anchoring model in Eq. 1 captures the anisotropic nature of liquid crystal films through the surface unit normal k. In biological membranes, the compressive strain and tensile lead to a contribution that is dependent to the gradient of k. Such mechanical effects are found in cell membrane (Lamparter and Galic, 2020;Jia et al., 2022), human neurons (Bianchi et al., 2019), mitochondria (Bostwick et al., 2016) and others Bagnani et al., 2021;Harbola et al., 2021). The Helfrich elastic model was applied to study the shape of red blood cell (Martínez-Balbuena et al., 2021) and lipid bilayer structure (Campelo et al., 2014). The elastic contribution to the anisotropic surface tension (Rofouie et al., 2017a) yields a membrane-liquid crystal tension (Rey, 1999b;2003;Shams et al., 2015) where κ is the surface curvature and k c the bending rigidity (Helfrich, 1973). By introducing the compression T 0 γ 0 − (k c /2)κ 2 0 , Eq. 2 reduces to a second order ordinary differential equation. Figure 6A and b present the daisy Bellis perennis and the SEM micrograph. Smaller corrugations are observed on top of a larger wrinkling peak, unlike Figure 2A. This delicate wrinkling pattern implies that another wrinkling mechanism should also play a role besides pure anchoring model. Figure 6C to f show the numerical simulations by introducing the elastic contribution into our shape Equation 4. Rofouie revealed that the chiral liquid crystalmembrane model is similar to a driven pendulum (Rofouie et al., 2017a), and four critical parameters were introduced to analyze the multiscale wrinkling patterns shown in Figure 6A and b:  (2017), from (Schauer et al., 2017); permission conveyed through Copyright Clearance Center, Inc. (C-F) are reproduced from (Rofouie et al., 2017a) with permission from the Royal Society of Chemistry.

Frontiers in Soft Matter
frontiersin.org l chiral and l membrane are the two length scales. ω is the dimensionless winding number, and the dimensionless W is the ratio of anchoring effect to mechanical compression. The pair (ω, W) can be used to categorize different models: pure anchoring model (CLC) implies that ω → ∞ and W W/γ 0 ; pure elastic model M) implies that ω W 0; and anchoring-elastic model (CLC-M) has finite values for both ω and W. The limit cycles and phase plots (as used in a pendulum model) for the three models are demonstrated in Figure 6E and f. Clearly, the elastic-liquid crystal model creates a complex wrinkling profile due to the constructive interaction of stress load and anchoring.
The wettability of multiscale wrinkling patterns on an elastic membrane is commonly studied (Jun et al., 2017;Li et al., 2019;Liu et al., 2022) due to its influence on surface chemical and physical properties (Vatansever et al., 2012;Deng et al., 2019;Zhou et al., 2020). Figure 7B and c demonstrate how the contact angle of the water droplet varies with different surfaces. The contact angle increases by 20.6°from an unwrinkled surface to a wrinkled surface, which validates that the contact angle is a tunable parameter that is dependent on the wrinkling profile. A significant control parameter is the dimensionless ratio of the height over spacing in the buckles (H b / S b ) (Byeon et al., 2020). reveals a linear relationship between the contact angle and H b /S b . The contact angle increases up to~60°with varying H b /S b . In partial summary, desired surface wettability can be achieved by adjusting the surface materials. To study the complex 1D open surface wrinkling, biomimetic approaches such as beetlesderived polarizers, and fish-skin inspired swimming suits can be widely applied in environmental-friendly industries.

2D open surface wrinkling: Egg carton corrugations
The egg carton surface is a ubiquitous nature's and engineering pattern that appear in various materials and processes, such as the electropolished zirconum (nm-scale), the papillae of rose petals ( Figure 8A) and dactyl club (µm-scale), erosion-driven rocks (cmscale) and others (Abukhdeir and Vlachos, 2011;Seeber et al., 2011).
Eq. 4 is the governing equation to the anchoring-driven wrinkling phenomenon in a 2D open surface. By applying a periodic boundary condition, Eq. (4) leads to a highly non-linear second-order partial differential equation. Using a perturbation method by expanding the surface energy and unit normal as follows  γ γ 0 + ϵγ 1 ( ) + . . . , k δ z + ϵk 1 ( ) + . . . where ϵ is a small parameter, δ z is a unit vector perpendicular to the (x, y)-plane and upper index 1) represents the first-order correction which leads to a linear shape equation. The solution to the shape equation reduces to a novel and revealing linear equation that connect geometry to orientation: where H is the mean curvature. Eq. 9 contains significant information: anchoring-driven wrinkling profiles along different directions can be superposed to create a complex wrinkling pattern. In Section 2, uniaxial wrinkling along direction z is simply sin qz. An egg carton surface naturally arises if two wrinkling profiles along orthogonal direction (x + y) and (x − y) superpose, as shown in Figure 8C and d, which verifies the egg carton structure observed on papillae of rose petal in Figure 8A and b. A 2D surface egg carton wrinkling is the result of balancing the dilation pressure, rotation pressure and director pressure as shown in Eq. 4, which shares the same wrinkling mechanism as a 1D surface discussed in Section 2.1. The second order derivatives of wrinkling profiles contains information of principal curvatures κ 1 < κ 2 . The mean curvature H = (κ 1 + κ 2 )/2, deviatoric curvature D = (κ 2 − κ 1 )/2 > 0, Gaussian curvature K = κ 1 κ 2 , Casorati curvature C (κ 2 1 + κ 2 2 )/2 and the normalized, dimensionless shape parameter S = (2/π) arctan (H/D) are the primary parameters to categorize basic shape information (Koenderink and Van Doorn, 1992;Wang et al., 2020b). By integrating information from these five measures, a complete characterization is made beyond the classical one based only on (H, K). Figure 9A shows the curvature lines for three primary shapes: sphere (magenta), cylinder (cyan) and saddle (yellow). An egg carton surface contains rich curvature features that can smoothly vary between those three primary shapes. Figure 9B is solved from Eq. 4 without linear approximation. The solution surface is composed by 5/ 31 egg carton, 10/31 cylinder, and 16/31 asymmetric biaxial wrinkling contribution. The non-linearity provides a more complex curvature distribution, shown in Figure 9B.
The control of local shapes has a prevailing application in cell biology (Malheiro et al., 2016;Werner et al., 2017;Bade et al., 2018;Pieuchot et al., 2018;Werner et al., 2018;Callens et al., 2020). Figure 9C and d are the schematics and experimental results of how geometry affects cell and tissue organization.
For a concave surface, the cell tends to be stretched, showing an upward stretched cell morphology. However, a convex surface allows cells to fully spread all over the surface, as shown in Figure 9C. A local cylinder has one vanishing principal curvature, providing two orthogonal directions for cell organization. Fibroblasts and mesenchymal stromal cells tend to avoid curvature and orient along the long axis (Hwang et al., 2009;Soiné et al., 2016;Bade et al., 2017;Werner et al., 2019), while epithelial cells prefer to wrap around the cylinder (Rovensky and Samoilov, 1994;Svitkina et al., 1995;Yevick et al., 2015;Werner et al., 2018), shown in Figure 9D ( Bade et al., 2018). discussed the cell organization on a saddle surface,  (Feng et al., 2008). Copyright ⓒ 2008 American Chemical Society. (C) and (D) are reprinted from , Copyright (2021), with permission from Elsevier.

Frontiers in Soft Matter
frontiersin.org creating more complex scenarios where no clear conclusions have been made (Bade et al., 2018;Werner et al., 2019;Callens et al., 2020). Cell migration, rather than static alignment, has also been observed to be highly dependent on the curvature of the substrate. For a convex cylinder (H > 0), hBMSCs migrate increasingly for the increasing curvature, while those on a concave cylinder (H < 0) show non-aligned mode and no angular preference Werner et al. (2018;2019). For a spherical surface, fibroblasts and mesenchymal stromal cells (MSCs) migrate faster if H < 0 Werner et al., 2017). Figure 9E presents the result of curvature-induced cell migration on an egg carton surface. The arrows show the tendency of cells to move with peripheral protrusions while maintaining a close distance of their nuclei to surface minima, see Figure 9F (Pieuchot et al., 2018). Figure 9G shows that the cell's nuclei are closer to the nearest surface minimum than the cell center, and the direction of the cell-nucleus vector always leads to the surface minima. The nuclei of cells tend to avoid the convex region and concentrate more on the concave region, exhibiting three times higher densities than flat region (see Figure 9H). The nucleus is a viscoelastic solid stiffer than the cytoskeleton (Guilak et al., 2000;Dahl and Kalinowski, 2011;Lammerding, 2011;Simon and Wilson, 2011). Experiments show force transmission between nucleus and cytoskeleton (Crisp et al., 2006;Dahl et al., 2008). For example, the deformation of the nucleus affects chromatin and further alters gene expression (Thomas et al., 2002;Cho et al., 2017;Szczesny and Mauck, 2017). Higher osteocalcin levels in MSCs were observed on a spherical surface, while on a planar surface, higher cytoskeletal forces were generated for osteogenic commitment (McBeath et al., 2004;Lammerding et al., 2006;Callens et al., 2020). The alignment and migration of cells are sensitive to the local curvature distribution (from Figure 9), therefore, cell nuclei are ideal curvature sensors (Wang et al., 2009;Anselme et al., 2018), and the control of curvature allows the engineering design for cell migration (Lauffenburger and Horwitz, 1996;Danuser et al., 2013;Yevick et al., 2015;He and Jiang, 2017; ), Copyright (2021, with permission from Elsevier. (B) Is reprinted from (Wang et al., 2022). Copyright (2020) by the American Physical Society. (C) and (D) are reproduced under CC-BY-4.0 (Callens et al., 2020). (E-H) are reproduced under CC-BY-4.0 (Pieuchot et al., 2018).

Frontiers in Soft Matter
frontiersin.org et al., 2017) and differentiation (McNamara et al., 2010;Eyckmans and Chen, 2014;Marino et al., 2014;Kim et al., 2015;Dobbenga et al., 2016;Lo et al., 2016). In partial summary, two orthogonal wave vector solutions to the liquid crystal shape equation generates uniaxial, equibiaxial and biaxial egg cartons with a rich spatial distribution of (H, K, C) curvatures and shapes S) that can be tuned and leveraged for biological applications such as tissue engineering due to the potential of the outstanding control of cell growth, differentiation and motion.

1D closed surface wrinkling: Droplets and tactoids
In this section, we present a review on the geometry of liquid crystal droplet enclosed by an elastic membrane. This structure is commonly discovered in bilayer lipid vesicles (Lawaczeck et al., 1976;Mueller et al., 1983;Liu et al., 2020), monolayer vesicles (Kunitake et al., 1981;Bauduin et al., 2011;Nam et al., 2019), chiral rafts (Sharma P. et al., 2014) and lyotropic liquid crystals (Iwashita and Tanaka, 2006;Tortora and Lavrentovich, 2011). (Guevara et al., 2019) studied four types of lipid-base droplets. Lipid-based droplets are essential carriers for mRNA encapsulation, which plays an important role in mRNA delivery. Lipoplexes are formed by cationic liposeomes and the phosphate backbone of mRNA (Elouahabi and Ruysschaert, 2005;Wasungu and Hoekstra, 2006;Guevara et al., 2020). Lipoplexes sometimes perform high instability and low transfection efficiency, and lipid nanoparticles are promising alternatives (Guevara et al., 2020;2019). The lipid nanoparticle consists of pH-responsive or cationic lipids enclosing polyanionic mRNA (Guevara et al., 2020). Lipid-polymer hybrid nanoparticles can be both bilayer and monolayer, showing different functionalities (Islam et al., 2018). The hybrid mRNA delivery organization has a higher nuclei acid condensation efficiency and tunable surface with multifunctionalities (Guevara et al., 2019). Figure 10A shows a phospholipid-coated liquid crystal droplet. Compared with lipid-based droplets (Guevara et al., 2019;, the surfactant-coated liquid crystal droplets have a larger scale (around 20 µm). The structure in Figure 10A can be used for biosensors, such as detecting antigens, pathogens, proteins, and antimicrobial peptides (Lowe et al., 2010;Bao et al., 2019). The bulk-surfactant interaction affects the surface orientation of the director (Brake et al., 2003a;. Especially, it has been shown that for such droplet, a critical length R~K/W determines the dominating factor between elasticity and surface anchoring, where surface anchoring energy~WR 2 and elastic energỹ KR He et al. (2014); Gupta et al. (2009);Miller and Abbott (2013). Besides the bulk-interface interaction, the surface of surfactant-coated liquid crystal droplet also responds to the environment such as pH, which further controls the director orientation inside the bulk (see Figure 10B. Rofouie et al. follow the theory from (Rey, 2000b;2006b;a, 2004b;a, Rey, 2005;Gurevich et al., 2014). Introducing Eq. 6 into Eq. 2 and assuming a non-vanishing pressure difference p c between the bulk and outer surface, the numerical solution solved with a periodic boundary condition results various morphologies of 1D droplet with an elastic membrane. The numerical solutions are given by Figure 11. Two dimensionless parameters are used in the simulation: the scaled anchoring coefficient 0 < |W/γ 0 | < 2 and the bending elasticity number 0.05 < α p c R 3 0 /k c < 50 (Meister et al., 1996;Ghochani et al., 2010;Boal, 2012), where R 0 is the spontaneous radius and is a constant R 0 = 0.5 μm in the simulation (Rofouie et al., 2017b). The pressure bending number β γ 0 R 2 0 /k c measures the ratio of the surface tension to the bending elasticity, while anchoring bending number ω WR 2 0 /k c measures the ratio of the anchoring to the bending elasticity. Figure 11 presents the membrane shape with (W ≠ 0) or without (W = 0) anisotropy. 2-fold (a)), 3-fold (b)), 4-fold (c)) and 5-fold (d)) symmetry are observed by varying the bending elasticity number α and the pressure bending number β. The morphology symmetry breaks when the anisotropic anchoring (W ≠ 0) is introduced. The curvature on the top of the membrane increases when the anchoring is enhanced. Comparing a) to d), the anchoring effect becomes less important in the symmetry breaking if the morphology shows more folds. Figure 11A is the stomatocyte shape observed in red blood cells (Tachev et al., 2004;Hajjawi, 2013;Muñoz et al., 2014;Geekiyanage et al., 2019). Figure 11E presents the relationship between the scaled bending energy and the reduced area for α = 0.3, 0.5, 0.7 and 1. A larger anchoring effect (higher W/γ 0 ) results in more folds, which enhances the scaled bending energy due to higher curvature. Figure 11f) demonstrates the effect of bending elasticity. A decreasing bending elasticity results in decreasing folds, where the morphology undergoes transition from 5-fold starfish to discocyte, ellipsoid, stomatocyte, umbonate (Chen et al., 2005;May et al., 2007), umbilicate (Sakamoto et al., 2017) and undulate (Suhandono et al., 2016), which are commonly used to describe red blood cells and bacterial colonies. A vanishing k c yields a spindle-like shape. The transition follows the curve in Figure 11F. In the end, by combining the three surface effects: bending elasticity k c , surface tension γ 0 and anchoring coefficient W, the morphologies are summarized in Figure 11G. Increasing k c and reducing γ 0 raises the number of folds while increasing W intensifies the asymmetry breaking on the top and bottom of a droplet.
In partial summary, simple extensions of the liquid crystal shape equation that include membrane elasticity generate morphology phase diagrams that include most of the observed cell shapes beyond the usual liquid crystal tactoids (Bagnani et al., 2018;Nyström et al., 2018;Gârlea et al., 2019;Li et al., 2020b;   Almohammadi et al., 2022). Some dynamics simulation of selfassembly can be found in (Khayyatzadeh et al., 2015;Han et al., 2017;Fu and Abukhdeir, 2018).

Conclusion and outlook
We presented a systematic review of recent experimental and theoretical literature of surface wrinkling patterns and morphology, and bioinspired applications of liquid crystal analogues. Table 1 summarizes the key aspects of different physical models (second row), wrinkling mechanisms and material parameters (third and fourth rows) and functionality (fifth row) of various biological systems that were discussed. The table shows the diversity of surface patterns and observed functionalities that can be captured by anisotropic soft matter capillarity.
Tulip Queen of the Night and spotted tinselfish show uniaxial surface corrugation, governed by the anchoring effect due to the preferred orientation of the cellulose and chitin fibres. The roughness of tulip petals and spotted tinselfish determines their optical and friction behaviour. The iridescence and glossy colour attract insects and further benefit pollination. In addition, an optimized surface roughness favours better swimming dynamics. The humidity of the environment changes the helix pitch of beetles' exoskeleton and vascular plants. The varying helix pitch disturbs the structural colour and diffraction grating, which makes beetles natural humidity detectors. Due to the significant role played by the elastic interface, the daisy flower petals show more complex multiscale corrugations. Both the anchoring and elasticity of petals govern the wrinkling patterns of daisy and kalanchoe blossfeldiana. Experiments validate that such a complex wrinkling profile enlarges contact angle, and it shed light on bioinspired hydrophobic materials. In partial summary, elastic biomimetic membranes governed by surface anchoring show complex wrinkling patterns, where the anchoring strength, elastic modulus and helix pitch are the dominant factors that control self-assembly and surface geometry.
Rose petals and dactyl club exhibit egg carton surfaces, where the non-vanishing Gaussian curvature contains rich geometric information. The egg carton surface of the biocompatible scaffold has been widely studied due to its promising application in controlling cell migration and cell differentiation. In the end, we also present a discussion on the morphology of lipid vesicles and surfactant-coated liquid crystals. A droplet displays complex anisotropy and morphology due to topological constraints. The curvature information of the elastic interface supports various biological applications such as mRNA encapsulation and virus/protein sensors. In partial summary, a tactoid with an elastic membrane is subjected to the topological constraint and performs an asymmetric profile that is highly dependent on the bulk-surface correlations.
In conclusion, liquid crystal capillary science, self-assembly and self-organization, and surface pattern formation were integrated to replicate, characterize, and shed light on many biological surfaces and functionalities. The outlook of this emerging field is rich in opportunities, which include surfactants, flow and mass transfer shaping, and reaction-polymerization-aggregation surface morphing, bioinspired diffraction gratings and structural colours, tunable and biodegradable scaffold for cell differentiation and tissue growth.

Author contributions
ZW, PS, and AR contributed to conception of the study. ZW wrote the first draft, prepared the figures and bibliography. ZW, PS and AR revised the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding
This work is supported by the Natural Science and Engineering Research Council of Canada (NSERC) (Grant number: 223086). AR is thankful to McGill University for financial support through the James McGill Professor appointment. ZW is thankful to the Faculty of Engineering for MEDA scholarship program to assist his PhD study.