Theoretical Platform for Liquid-Crystalline Self-Assembly of Collagen-Based Biomaterials

The collagen triple helix is a ubiquitous biomacromolecule used in acidic aqueous solutions as precursor in the fabrication of artificial compact bone and cornea and in tissue engineering. The primary architecture of these highly structured solid tissues is formed during the cholesteric liquid-crystalline stage of their morphogenesis. The theoretical platform that describes the coupled dynamics of phase-ordering and mass transfer developed, implemented and validated here can be used for optimal material design and plays a significant complementary role to future experimental studies. Based on uniaxiality assumption, we have recently developed and validated a theory for the free energy tailored for acidic collagenous dispersions. Here we significantly expand and generalize our previous study, by including biaxiality since cholesteric phases must have a degree of biaxiality. In this work, we first modify the proposed interchain interaction and excluded-volume contribution by use of the addition theorem for spherical harmonics. Then, the Euler-Lagrange minimization followed by expansion around I/N* transition allows us to construct the free energy of ordering in terms of the phenomenological Landau–de Gennes formulation. Finally, we use the time-dependent Ginzburg-Landau equations to study the non-Fickian evolution of a single two dimensional cholesteric tactoid through a shallow quench from the isotropic to biphasic region of the phase diagram. Although equilibrium biaxiality is considerably low for these long-pitch cholesterics, we found that during self-assembly the biaxial order parameter achieves significant larger values than the equilibrium value. Additionally, the relaxed director field becomes both onion-like and defect-less, which is consistent with the twisted bipolar structure observed experimentally. The self-assembly simulations demonstrate that the formulated theoretical platform is not only consistent with previous theoretical and experimental studies but also able to be used to explore new routes for non-equilibrium collagen self-assembly. Taken together, this study deepens our understanding of cholesteric (chiral nematic N*) mesophase in acidic solutions of tropocollagen, and suggests a systematic spatio-temporal model that is capable of being used to extract the engineering principles for processing of these sought-after biomaterials.


INTRODUCTION
Type I Collagen is composed of three left-handed polypeptide helices (denoted by [α1 (I)] 2 [α2 (I)]) twisted together to yield a right-handed triple helix. This rod-shaped biomacromolecule, also known as tropocollagen, commonly has a 1.5 nm bare diameter and 300 nm height. The tropocollagen falls into the class of fibrous proteins and is abundantly found in both soft and hard human's tissues, namely cornea, tendon, cortical bone, and more [1]. Over the past two decades, biomimetic fabrication of collagen-based biomaterials has received considerable attention in view of the abundant critical applications such as artificial bone [2][3][4][5] and cornea [6,7] reconstruction. Moreover, for invitro replication of these collagenous biological tissues, there is fortunately no concern about supply because tropocollagen can be readily accessible through mammalian and non-mammalian resources [8]. Consequently, numerous promising applications of biomimetic fabrication of collagenous biomaterials [9][10][11] in conjunction with the availability of precursor play a central role in the drive to create the bioinspired collagen-based materials.
The structural pattern of tropocollagen rods bestows great structural-relation properties on collagenous biological materials and biomaterials. Furthermore, their structures are analogous with architecture of tropocollagen in liquid-crystalline states [12], hence these materials are called "solid analogs." This correspondence establishes the role and impact of liquidcrystalline morphogenesis [13][14][15] and singles out liquid-crystalbased biomimetic material process engineering as a promising route to enhance the quality of collagen-based biomaterials or even to explore new ones [3,[16][17][18][19].
Normally, tropocollagen is immiscible in aqueous solutions due to its hydrophobicity. To attain a stable aqueous isotropic phase, which is the starting point of biomimetic fabrication, hydrophobicity of tropocollagen must be reduced by being dispersed in acidic solutions. Basically, numerous amine function groups that are good proton receptors are found along the tropocollagen backbone. Once these functional groups are protonated, the intrachain repulsion causes that, the semi-flexible (worm-like) backbones become uncoiled and essentially rigid rods. The existing interchain repulsion also impedes aggregation, in other words the rods have an effective diameter between two or three times the bare one [20][21][22]. Finally, due to being chargecarrier rigid rod-like molecules, tropocollagen is capable of exhibiting lyotropic cholesteric phase organization. For example, for an acetic acid concentration of [AC] ≈ 2, 900 mM, a phase transition from isotropic to chiral nematic (N * ) takes places at tropocollagen concentrations of [C] ≈ 88 mg/ml [21].
Although the primary architecture of these versatile biomaterials are formed at the molecular level (i.e., mesophasic stage), the focus has been at the tissue level [23] and studies on molecular level are few [21,22]. Furthermore, to the best of authors' knowledge, theoretic studies of cholesteric selfassembly of aqueous acidic tropocollagen solutions have not been carried out, which also reflects the case of chiral nematic phase ordering in general [24][25][26]. To address this gap, we have recently developed, implemented, and validated a theoretical model tailored for self-assembly of tropocollagen dispersed in acidic aqueous solutions [20]. This thermodynamic theory [20], which is based on the uniaxiality assumption, has integrated microscopic mechanisms of mixing entropy and enthalpy, attraction, repulsion, twisting, excluded-volume, and chirality. In the present study, we lift the uniaxiality assumption by generalizing the free energy that includes biaxial effects. This is crucial for cholesteric materials because chiral nematic phase is described by two vectors: the director (n) and the chiral axis (h), additionally cholesterogens are intrinsically biaxial as discussed by Wulf [27] and Wright and Mermin [28]. Incorporation of the biaxial order parameter into the cholesteric self-assembly deserves consideration because biaxiality influences patternformation even in nematic mesophase, such as interfacial biaxiality under tangential director orientation [29][30][31][32], the biaxial core of singular disclinations [31,33], and sometimes more pronounced biaxiality under time-dependent conditions than under static equilibrium [34]. For the above reasons we first include biaxiality in the model formulation stage and then focus on its emergence in bulk, defect core, and interfacial regions; which are of significant importance in all structured materials [30,33,[35][36][37].
In our previous validated work [20], we showed that our thermodynamic model of acidic collagen solutions captures two key features: (i) the expected chimney diagram predicted by Flory and found experimentally for many lyotropic rod-like liquid-crystalline polymers [38], and (ii) the parabolic bi-phasic funnel in aqueous acidic collagenous solutions under increasing pH, where cholesteric tactoids (drops) emerge from isotropic phases. Study of cholesteric tactoids is important because of three main reasons: (i) tactoid formation process must occur in to chimney and funnel phase diagrams, which are the fingerprint of rod-like macromolecules. Thus, these cholesteric drops are a crucial element in the validation of thermodynamics of rodshaped rigid macromolecules; (ii) these stable but deformable drops serve a sources of material properties information such as bulk Frank-Oseen-Mermin elasticity [27], novel coupled gradient contributions between nematic order parameter and collagen concentration, and the cholesteric pitch; (iii) characterizing and understanding the emergence, growth, annihilation, and coalescence of tactoids are essential to future developments of collagen-based material processing. To focus on collagen tactoids, as shown in Figure 1, we then target the dynamic of selfassembly through a shallow quench from isotropic phase into the bi-phasic funnel of the previously obtained phase diagram [20]. In contrast to the better known single component monomeric thermotropic tactoids, in the present case concentration is a conserved transport variable that need to be included. For this purpose, we formulate the coupled phase ordering/mass transfer Model C [39,40] in order to derive the governing equations of collagen self-assembly. Afterward, we impose proper auxiliary conditions (e.g., initial and boundary conditions for the computational domain) on the obtained governing equations to capture a thorough spatio-temporal evolution of a single cholesteric tactoid-see Figure 2. This evolution has two steps: (i) emergence of a cholesteric nucleus in a continuous isotropic phase, (ii) followed by the formation of a stable chiral nematic tactoid coexisting with the isotropic phase.  [20]. The schematics denote the isotropic phase at low collagen concentrations, a typical micron-sized cholesteric drop in an isotropic bulk at intermediate concentrations, and the chiral nematic (N*) or cholesteric phase at higher concentrations.
In this work we restrict simulations to a single collagen tactoid with the aim of contributing to the evolving understanding of chiral phase ordering [24,25,41]. The simulations are also restricted to 2D. In principle, 3D spatio-temporal simulations can give a full picture of tactoid formation stages. Yet, from practical viewpoint, the present phase ordering/mass transfer coupled non-linear model with nano-to-micron scales becomes essentially intractable [42]. Furthermore, we have previously shown [43][44][45] that 2D simulations can provide invaluable predictions, and as discussed later on, in this study the important metrics of size, shape, and structure are not lost when using our 2D simulation box. In particular, we capture bulk disclinations, interfacial anchoring, interfacial biaxiality, growth modes, and self-selected shapes. Hence, this 2D study gives a necessary foundation for future 3D simulations.
The paper is organized as follows. Section Continuum Methodology for Simulation of Liquid-Crystalline Self-assembly of Tropocollagen Dispersed in Acidic Aqueous Solutions presents the methodology used in the formulation of self-assembly, including: (1) Formulation of the free energy for a system consisting of charged cholesterogen dispersed in a mixture of water solvent and mobile ions-see subsections Long-Range Description of Molecular Alignment to Total Free Energy Tailored for Tropocollagen Self-assembly in Acidic Aqueous Solutions. Subsection Long-Range Description of Molecular Alignment defines the Q-tensor. In Free Energy Contributions for Pure Charged Cholesterogens; Incorporation of Biaxial Order Parameter, the free energy of pure charged cholesterogen is developed taking into account the biaxial order parameter. In Mixing Free Energy of Binary Dispersions Consisting a Charged Cholesterogen and Small-Sized Solvent, the obtained free energy is generalized for a mixture of charged cholesterogen and smallsized solvent. In Total Free Energy Tailored for Tropocollagen Self-assembly in Acidic Aqueous Solutions, we discuss and incorporate other free energy contributions involved in the evolution of mesophasic state, such as the elasticity of Frank-Oseen-Mermin [28] and gradient contributions, and formation of the I/N * interface. Thus, in Total Free Energy Tailored for Tropocollagen Self-assembly in Acidic Aqueous Solutions, the total free energy of system is formulated [2]. Formulation of governing equations along with the appropriate auxiliary conditions for simulation of liquid-crystalline self-assembly in which a cholesteric nucleus of tropocollagen is initially seeded and allowed to spontaneously growth in coexistence with isotropic phase-see subsections Governing Equations for Kinetics of Self-assembly; Orientational Relaxation, and Uphill Diffusion to Computational Details. In Governing Equations for Kinetics of Self-assembly; Orientational Relaxation, and Uphill Diffusion, the governing transport equations (Model C) are formulated. Finally, subsection Computational Details presents the implementation of self-assembly simulation for nucleation and growth of a single cholesteric tactoid coexisting with an isotropic phase. Section Results and Discussions, presents results of emergence and growth of a cholesteric tactoid. Lastly, the conclusions and nomenclature are summarized in sections Conclusions and Nomenclature, respectively.

Long-Range Description of Molecular Alignment
The long-range orientational order in a liquid-crystalline phase is parameterized by a second-order symmetric traceless tensor called Q-tensor [14,15,46].
where δ is the Kronecker delta. The orientation of a mesogen is characterized by the orthogonal director triad of (n, m, l). The degree of alignment along the uniaxial director, n, and biaxial director, m, are S and P, respectively. Due to the quadrupolar symmetry of Q-tensor, it possesses the salient feature of head-tail invariance of molecular alignment (i.e., n ≡ −n, m ≡ −m, and l ≡ −l). The largest absolute eigenvalue of Q-tensor equals to 2S/3 and the corresponding eigenvector is equivalent of uniaxial director, n. The difference between the absolute medium and smallest eigenvalues is 2P/3 and the eigenvector corresponds to the second largest absolute eigenvalue is biaxial director, m. Thus, in the isotropic and ordered phases, the Q-tensor becomes the 3 × 3 zero matrix, Q = 0, and non-zero matrix, Q = 0, respectively. The uniaxial and biaxial order parameters are also defined in terms of directors/Q-tensor or the normalized orientational distribution function on the unit sphere, ψ (u), for any given molecular orientation, u: dθdϕ represents a solid angle, and θ and ϕ are the polar and azimuthal angles. ( θ , ϕ) is defined as 3sin 2 (θ) cos (2ϕ)/2. As explained below, P 2 (cos (θ)) and (θ, ϕ) are representative of uniaxiality and biaxiality, respectively. In addition, the normalized distribution function, employed in Equations (2, 3), implies following constraint [47]:

Free Energy Contributions for Pure Charged Cholesterogens; Incorporation of Biaxial Order Parameter
The total dimensionless Helmholtz free energy per particle, F, for a dispersion comprising N A charged cholesterogens is Frontiers in Physics | www.frontiersin.org [20,48,49]: where β, µ o (T), and c A stand for thermal energy, standard chemical potential and number density, respectively. The last three terms in Equation (5) account for the contribution of molecular orientation (i.e., uniaxiality and biaxiality) in the mesophasic state. σ (ψ (u)) describes the decrease of orientational entropy upon alignment of the mesogenic molecules.
Since we focus on rod-like rigid mesogen, the second virial approximation is capable of accurately describing the excluded volume effect, given by Odijk [49] and Drwenski et al. [50]: υ AA is the average excluded volume defined as πD eff L 2 /4 in which L and D eff denote contour length and effective diameter of tropocollagen. D eff has a dependence on the bare diameter, D = 1.5 nm, and concentration [20]. To take biaxiality into consideration, in accordance with Drwenski [51,52] and Matsuyama and Crystals [53], we make use of the addition theorem for spherical harmonics to express the angle between two rods, γ, in terms of the polar, θ, and azimuthal, ϕ, angles in spherical coordinate: P 2 (cos (γ)) = P 2 (cos (θ)) P 2 cos θ ′ + (θ,ϕ) (θ ′ , ϕ ′ )/3 (8) First term in Equation (8) is independent from azimuthal angle and represents uniaxiality. Second term is related to biaxial contribution and has the dependence on both polar and azimuthal angles. The intermolecular interaction and angle between rods interchangeably affect each other because the electrostatic repulsion and twisting favors perpendicular orientation while the van der Waals attraction prefers the parallel alignment (i.e., nematic phase) [49,50,54]. Hence, based on our previous work [20] and [51], we suggest the net interchain potential expressed by where υ A is the volume of an individual rigid rod, υ A = πD 2 eff L/4. U ′ elc and U ′ MS are parameters of electrostatic repulsion and a positive constant, respectively. The contribution of intermolecular interaction, M (ψ (u)), is then obtained by taking average over all possible rod configurations [20]: U = U elc − U MS is called the potential of the orientationdependent intermolecular interactions where U elc and U MS are the strength of electrostatic repulsion and Maier-Saupe constant that is a positive constant independent of temperature. Note that Q : Q is related to uniaxial and biaxial order parameters by 2 S 2 +P 2 /3 /3. We note that the effective diameter reflects the intermolecular repulsion, or to put it another way, the effective thickness of the attached ions on the backbone of tropocollagen. This effective thickness is called double-layer thickness, ακ -1 [49,50]: α and κ -1 , which are defined as follows, are parameter of doublelayer thickness and Debye screening length, respectively: where λ B , N avo , ǫ, m, Z, Ei (•), γ E , and are the Bjerrum length, Avogadro's number, ionic strength, molar concentration, charge number, the exponential integral defined as Ei (x) = − ∞ -x exp (-t)/t dt, Euler constant equals to 0.5772, and linear charge density. A detailed account of parameters' values, their selection and physical significance and physical properties for aqueous acidic collagen I solutions is given in Khadem and Rey [20].

Mixing Free Energy of Binary Dispersions Consisting a Charged Cholesterogen and Small-Sized Solvent
The mixing free energy of the binary solution is given by Matsuyama and Kato [48] where F s (N A ,N I ), F s (N A , 0) and F s (0,N I ) are free energies of solution, pure anisotropic component dispersed in isotropic state and isotropic component, respectively. Thus, in this subsection, we shall first derive the free energy of solution, and then formulate the mixing free energy of a binary dispersion by use of Equation (13).
Substituting Equations (7, 10) into Equation (5) leads to the free energy of pure charged chiral nematic rods. The free energy for binary mixture of charged chiral mesogen and smallsized solvent (water in our case)-denoted by subscript A and I, respectively-is then formulated as Equation (14) is not usable unless the unknown normalized distribution function, ψ (u), is known. To formulate the normalized distribution function, the total free energy of system subjected to the normalizing constraint, given by Equation (4), is minimized using Euler-Lagrange method. This minimization yields an irreducible algebraic integral equation expressed by Simplicity of free energy expression is essential since our ultimate objective is the self-assembly simulation which in itself is computationally complex. A heavy computational load is expected because the self-assembly process covers a wide range of length scale (i.e., ranging from nano-to macro-scale) and it may go through a variety of complex microscopic mechanisms [41,[55][56][57][58][59]. Thus, to improve tractability, we expand the functional part of Equation (15), Ŵ (γ), in terms of the second Legendre polynomial by use of Equation (8): Having substituted Equation (16) into Equation (15), the normalized distribution function is obtained: where the I 00 is a definite integral defined as (18) W is known as the net cholesteric potential, which is similar to Khadem and Rey [20] and can be parameterized as φ A is the effective volume fraction and h = (κD eff ) −1 . α W is assumed to only be dependent on concentration of acid throughout the evolution-a reasonable assumption because α W is mainly affected by concentration of acid [20]. Next the mixing free energy, Equations (20a-c), is obtained by use of Equations (13,14,17). Detailed account of such algebraic derivation are given in Khadem and Rey [20] and Odijk [48]. Note that hereafter, for convenience, we use φ to represent the effective volume fraction of tropocollagen-it can be related to the concentration of tropocollagen in units of mg of tropocollagen pre ml of solution by C = φ/α c where α c is a unit conversion factor. The dimensionless mixing free energy density is: where n stands for number of segments on tropocollagen backbone.f iso andf h describe different physics; the former explains the phase separation which is the well-known Flory-Huggins equation and the latter controls the phase transition (i.e., homogenous contribution). In the absence of biaxiality, P = 0, the obtained mixing free energy, Equations (20a-c), is reduced to the validated free energy functional given in Khadem and Rey [20] which was validated with experimental data of tropocollagen and with previous theoretical studies. It is worth mentioning that with further assumptions the obtained free energy density leads to the formulation given in Matsuyama and Kato [48] as well as the well-established theory of Onsager-see ESI of Khadem and Rey [20] for further discussion. For numerical tractability, similar to Matsuyama [51] and Matsuyama and Kato [48], we make use of a Taylor expansion in vicinity of I/N * to expand Equation (20c) in a power series of order parameters, S i P j -the resulting polynomial is the phenomenological Landau-de Gennes (LdG) theory [60]: Although self-assembly simulations by use of the Equations (21ad) is appreciably more tractable than with Equation (20c), it should be noted that the used expansion may affect the accuracy of simulations in the cases of deep quenches. However, this study only focuses on the self-assembly of shallow quenches into biphasic region which is a narrow region around I/N * boundary, see Figure 1.
An order-disorder phase transition takes place if and only if W = α W φ = α W α c C > 5 to make the coefficient of second invariant of Q-tensor, a, negative. The derived LdG coefficients satisfy two general theoretical expectations; (1) the first-order phase transition (i.e., B = 0), and (2) two minima correspond to isotropic and ordered phases (i.e., A<0 and C>0) [46]. Furthermore, in the cholesteric phase, W can be about 10, and under such conditions the proposed LdG coefficients becomes similar to the well-established lyotropic LCP Doi's model, were b ≈ c [14,15,61].
Total Free Energy Tailored for Tropocollagen Self-Assembly in Acidic Aqueous Solutions In addition tof mxing which is capable of describing phase separation and an order-disorder phase transition, for constructing the total free energy of mesogenic solutions, the contributions of gradients should be taken to account [62][63][64]: ξ = a 3 L 1 β is the coherence length in which a 3 stands for the volume of each lattice unit and L i are elastic constants. ∇ = h 0 ∇ is dimensionless gradient in which h 0 denotes a macroscopic length scale and the spatial domain is scaled by where L φ is cost of interfacial formation and L φ−Q represents coupling constant. The total free energy as well as the evolution of chiral nematic phase for tropocollagen are mesoscopic because it retains both microscopic length scale, ξ, in a nanometer range and macroscopic length scale, h 0 , in the range of micrometers.

Governing Equations for Kinetics of Self-Assembly; Orientational Relaxation, and Uphill Diffusion
Simulations of pattern-formation in fibrous composites, including collagen-based tissues, were first carried out by De Luca and Rey [61,64,65]. Their approaches were based on diffusionless evolution of mesophase, and capable of predicting macroscopic architecture of these materials to a great extent. However, recent studies have revealed the imperative role of diffusion in accurately capturing the growth of order-disorder interface [43,45]. Hence, for the purpose of realistic self-assembly modeling, in this subsection, we formulate the spatio-temporal evolution of tropocollagen in which the Q-tensor augmented with a mass transfer equation.
The cholesteric micro-structures in collagenous biomaterials are formed through the liquid-crystalline self-assembly stage.
Two simultaneous mesoscopic mechanisms govern this thermodynamically driven assembly. First, mass transfer mechanism allows tropocollagen macromolecules to diffuse into cholesteric phase (i.e., tropocollagen-rich phase) from isotopic phase (i.e., tropocollagen-lean phase). The mentioned demixing is known as uphill or non-Fickian diffusion and reduces the total free energy of system. Second, orientational relaxation mechanism induces cholesteric architecture inside the formed high-concentration domain. To describe these two phenomena; two coupled fields are required. First, the conserved scalar field of concentration, C, or equivalently volume fraction, φ , governing the phase separation. Secondly, the non-conserved tensorial field of Q-tensor by which the orientation of tropocollagen biomacromolecules is primarily specified. The spatio-temporal evolution of {Q, φ } is found using the time-dependent Ginzburg-Landau (TDGL) formalism, also known as model C in Hohenberg and Halperin classification [39,40]. The dimensionless form of model C adjusted for self-assembly simulation reads [17, 61,64,66,67]: δF net /δQ represents functional derivative.t is dimensionless time defined ast = tM Q / a 3 β where t is time,M φ = M φ / M Q h 2 0 in which the mobilities of Q and φ are M Q and M φ , respectively. Additionally, the superscript [s] indicates that the functional derivative must be symmetric traceless in order to be consistent with the nature of Q-tensor-for any given second rank tensor where superscript t denotes transpose.
The system given in Equations (23a,b) is a set of six coupled non-linear PDEs. Equation (23a) accounts for the spatiotemporal evolutions of the orientational tensor order parameter. This equation is the compact tensorial form of five independent second-order PDEs. Furthermore, Equation (23b) is a fourthorder PDE, known as the Cahn-Hilliard equation, to describe the concentration field by which the chiral nematic and isotropic phases gradually evolve through the uphill diffusion mechanism.

Computational Details
Here we elaborate on the simulation of nucleation and growth of an isolated cholesteric tactoid in a continuous isotropic phase. This simulation consists of a diffusional phenomenon coupled with structural relaxation. The general schematic representation of this implementation is illustrated in Figure 2.
As above mentioned, the biomimetic formation of collagenbased tissues starts with dissolving tropocollagen in acidic aqueous solutions to obtain the isotropic phase. In such condition, a nucleus is thermodynamically allowed to grow, providing its radius is greater than a critical value. In that case, as a single tactoid grows, the tropocollagen rods diffuse from collagen-lean phase to collagen-rich phase, in turn, the isotopic and cholesteric phases become depleted from and enriched in tropocollagen, respectively. The diffusion of tropocollagen from lean phase (isotropic phase) to rich phase (cholesteric phase) continues till a point where the chemical potentials of two phases become identical.
As illustrated in the Figure 2, we consider the bulk of system as a square with [−0.5 0.5] × [−0.5 0.5] normalized by h 0 . Each pair of sides are subjected to the periodic boundary condition. Initially Q = 0 and the phase is isotropic. Afterwards, an initial cholesteric tactoid is seeded by a circular Gaussian distribution with FWHM (full width at half maximum of Gaussian function) greater than the critical drop diameter. The seeding is expressed: and , that are, respectively, a second rank symmetric traceless random tensor and scalar random number, are included in the modeling to represent the fluctuations existing in a real system. The subscript e indicates the equilibrium condition given by: In accordance with l = n × m, l e is computed as 0 0 1 . Equations (24a,b, 25a,b) describe a nucleus whose center placed at position (x 0 , y 0 ) at a concentration equivalent to the effective volume fraction of φ ch . For convenience, we choose the center at (x = 0, y = 0). The concentration of tropocollagen from the center, which is a cholesteric phase, gradually decreases along the radius to the concentration of continuous isotropic phase, φ iso . This approach for simulating the initial nucleus was adapted from Wincure and Rey [30,31,33]. In order to make sure that the initial drop is sufficiently large, we choose the σ x = σ y = σ and obtain σ in way that FWHM equals two times the critical diameter: FWHM = 2 √ 2 ln 2σ = 2D c = 4R c . Classical Nucleation theory [68] provides a rough estimation of the critical radius of a drop as expressed by in which γ i and µ iso−Cho are the interfacial tension and the chemical potential difference between isotropic and cholesteric phase [68,69]. Additionally, Equation (24b) yields the quenched concentration as φ q = CD φ| t=0 dxdy/ CD dxdy in which CD denotes the entire system (computational domain). Consequently, φ iso plays an appreciable role in the size of tactoid because its value affects the initial amount of tropocollagen existing in the system. Furthermore, the total conservation of mass is imposed by: Simulation parameters used in this study are summarized in Table 1-also readers are referred to the Khadem and Rey [20] for detailed account of parameters selection in order to accurately capture the available experimental data. Although theL φ−Q ,L φ , and L 2 /L 1 have not been documented for tropocollagen, we choose common values which satisfy the energy transformation constraint [70][71][72]: Equations (24a,b) in conjunction with the above-explained conditions are solved with an adaptive finite elements technique with biquadratic basis functions (General PDE solver of COMSOL Multi-physics 5.3a). Furthermore, to acquire the acceptable spatial resolution, we considered at least 50 elements per pitch which resulted in nearly 10 4 triangular elements, and temporal resolution was carried out by the Backward Euler method. Convergence, accuracy, and stability were checked using standard techniques-for further information on the method and solution approach, please see the accompanying Supplementary Material.

RESULTS AND DISCUSSIONS
In this section, the dynamics of mesophasic evolution and the resulting equilibrium configuration for a shallow quench from the isotropic phase into the cholesteric phase in the presence of one small cholesteric seed are given and discussed (see Figures 1, 2). Figures 3a-d show snapshots of a growing tactoid corresponding to dimensionless times 0, 900, 960, and 1,100, respectively. Figure 3a shows the initial condition of a small chiral nematic drop seeded in a large isotropic phase area. Note that only a small section of the computational domain, in which the self-assembly is supposed to take place, is shown in Figure 3. The computational domain is actually chosen as a fairly large square with length of h 0 = 100 µm in order to make sure that the existing amount of tropocollagen in the system is sufficient for formation of a single cholesteric tactoid with diameter of the order of 30 µm-as experimentally observed [21]. The size of initial seed must be greater than a critical value in order for the drop to grow based on the mechanism of uphill diffusion, otherwise downhill diffusion takes place and the initial drop is dissolved in isotropic phase.
Although the initial configuration of rods is chosen as twisting around x-axis, see Figure 3a1, the rods prefer to be aligned in a concentric configuration, as shown in Figure 3b1. During the early growth of the tactoid, rods attempt to radially twistthe helicoidal axes are along the radii of the circular tactoid. Yet, rods placed at the center of drop exhibit orientational frustration. This frustration emerges int = 900, Figure 3b1, and yields a τ +1 cholesteric defect. As the tactoid grows, the central rods resolve the orientational frustration with an escaped configuration (see Figure 3d1) known as a non-singular λ +1 cholesteric disclination. These important predictions may be The square-brackets next to the values indicate the corresponding unit, and [-] shows dimensionless. For those parameters which have not been documented for solutions of tropocollagen, the common values are used instead. Readers are referred to Khadem and Rey [20], Gobeaux et al. [21], De Luca and Rey [61], Gurevich et al. [62], Das and Rey [63], and De Luca and Rey [64] for details of parameter selection.

FIGURE 3 |
The spatial distributions of order parameters, S and P, in conjunction with the director configuration at the early growth of cholesteric tactoid shown in (a-d). In the first column, the uniaxial configuration, n, of tropocollagen macromolecules are represented by rods whose color (blue to red) shows the uniaxial order parameter, S. To complete the understanding about the configuration of rods in xy-plane, in the second column, the z component of n is shown by use of a monochromatic blue spectrum. In last column, the monochromatic cyan denotes the variation of biaxial order parameter during the time evolution. (e) Illustrating the color bars for S, n z, and P, the used coordination of system and length-scale bar.
difficult to be captured experimentally due to intrinsic size length scale resolutions when using optical methods [73][74][75]. Figures 3a2-d2 show the z component of uniaxial director, n. The figures show that the central director regions evolves slower and lags the radial helix formation that results in tangential interfacial orientation experimentally observed for tropocollagen tactoids [37]. The tangential orientation minimizes the interfacial free energy at n.k = 0 where k is the interfacial normal vector. This tangential configuration, n.k = 0, emerges when the coupling coefficient,L φ−Q <0 [21,71,72,76]. The structure of the 2D tactoid is a radial helix, with tangential interface orientation at the edge and non-singular escape orientation at its center.
Of particular interest to this study is incorporation and analysis of biaxial order parameter during the evolution of the cholesteric tactoid. In the third column of Figures 3a3-d3, the spatial variation of biaxial order parameter, P, is shown in the early stages of growth. The biaxial order parameter becomes particularly noticeable at the interface and at the defect core. Thus, we found that although the equilibrium biaxiality for rod-like macromolecules is small [27,51], during the phase ordering it takes a larger value than its equilibrium; the difference between dynamical and equilibrium values for biaxiality may be up to three orders of magnitude.
In the course of time, the Q-tensor is relaxed, mass transfer ceases and the structure equilibrates, as shown in Figure 4. As depicted in Figures 4a1,b, the equilibrium configuration of tropocollagen rods becomes concentric, also known as onionlike. This defectless configuration, which has a non-singular λ +1 cholesteric disclination at its center, thoroughly matches with the xy-cross-section of Twisted Bipolar Structure (TBS) given in Sec et al. [77]. Moreover, the experimental POM image reported in Gobeaux et al. [21] confirms TBS for the 3D tropocollagen tactoids. Consequently, the resulting 2D configuration, shown in Figures 3, 4, is consistent with experimental observation. Figure 4a1 shows that the size of tactoid is also consistent with experimental results given in Gobeaux et al. [21]. As can be seen, the diameter of tactoid contains three pitches that each of which has a length of 10 µm. Therefore, the tactoid shape becomes a nearly 30 µm spherulite. Figure 4c represents the equilibrium concentration field. Although a gradient of concentration exists in the interface, the drop remains intact and stable in the isotropic phase. This feature verifies that the growth of cholesteric tactoid is according to the mechanism of uphill diffusion. Figure 4d demonstrates that the equilibrium biaxial order parameter P in the interface is nearly 0.04, however its value sharply decreases to 10 −4 confirmed by previous theoretical studies [27,51].
Through the entire evolution we find: (1) the interfacial uniaxial order parameter is approximately S c = 0.39 that is quite close to the critical uniaxial order parameter reported in Khadem and Rey [20] and Gobeaux et al. [21]; (2) The biaxial order parameter at the tactoid's interface is at all times greater than in the interior. The only exception is during initial defect nucleation where biaxiality pronouncedly appears at the core and edge of the 2D drop.
According to the principle of minimum free energy, the kinetic of a spontaneous process follows a path over which the free energy progressively decreases and ends up in a minimum at equilibrium. Figure 5 illustrates the averaged free energy contributions through the phase ordering process of tropocollagen dispersed in the constant concentration of 2.9 M acetic acid. These spatial averages are defined as The formation of the single cholesteric drop is the interplay of five free energy contributions. The entropic and enthalpic contributing factors in isotropic phase separation are described by Flory-Huggins theory,f iso . The LdG theory,f h , also accounts for the homogeneous effect of phase transition. The spatial averages of these contributions are shown bȳ F iso andF h , respectively. The monotonic decrease inF iso andF h shows that the phase separation and phase ordering are energetically favorable. In addition, it emphasizes on the lyotropic nature of phase ordering in acidic collagenous dispersions; rods are spontaneously accumulated in cholesteric phase, in turn, removed from the isotropic phase. Hence, F net which is the summation of all contributions, is considerably affected by contributions of phase separation and phase ordering. In spite of these energetically favorable contributions, formation of I/N * interface and cholesteric configuration inside the tactoid require energy costs which are reflected as penalty terms in the net free energy; see Equations (22a-d). The green solid line in Figure 5,F cg , depicts the cost of interface formation (i.e., mass gradient zone shown in the Figure 4c). This cost is nearly 40 percent of the energy reduction in either phase separation,F iso , or phase ordering,F h , thus the interfacial formation cost can be compensated. Furthermore, the black dash line,F e , and the purple dash line,F c , stand for the average costs for the onion-like configuration of rods inside the chiral nematic tactoid and the tangential configuration in interface, respectively. As seen, the formation cost of the interfacial parallel anchoring,F c , is ∼2% of the interior cholesteric energy,F e .

CONCLUSIONS
Building on our prior work [20], in this study, we have developed and validated a theoretical framework to study the spatiotemporal phase ordering of tropocollagen dispersed in acidic aqueous solutions into 2D drops. By use of the addition theorem for spherical harmonics (Equation 8), we first incorporated the biaxial order parameter P (Equation 3) into the orientational entropy (Equation 6), the second virial approximation Equations 7a-c), and the intermolecular interaction (Equation 10). We then obtained the LdG coefficients (Equations 21a-d), and formulated the net free energy of system, (Equations 21a-d).
To capture the kinetic of the emerging 2D tactoids size, shape, and structure, we relied on the proposed net free energy and phase ordering/mass transfer process (Model C) to establish the governing equations, Equations (23a,b), which were numerically solved under the mentioned auxiliary conditions elaborated in subsection Computational Details.
Figures 3a-d reveal two findings. First, the physical origin for the non-singular escaped λ +1 disclination. Basically, in the early evolution a τ +1 defect emerges at center of nucleus. As time progresses, the central directors go through a defect shedding stage and the τ +1 cholesteric defect evolves into the escaped λ +1 disclination. Second, at the interface and defect core region, the biaxial order parameter takes appreciably large value in the early evolution. Furthermore, Figures 4a1-a3,b demonstrate that the resulting equilibrium state of collagen tactoid is an ∼30 µm spherulite in which the rodshaped macromolecules are aligned in concentric configuration, consistent with experimental observations [21]. Taken together, these results contribute to the development of optimized processing protocols for collagen-based materials and their material property characterization.

DATA AVAILABILITY
All datasets generated for this study are included in the manuscript and the Supplementary Files.

AUTHOR CONTRIBUTIONS
SK developed the theoretical framework, carried out the simulations, and analyzed the results. AR supervised the study, analyzed and interpreted the results. All authors discussed the results, contributed to writing, and agree about the content.