C OMPLEX A RCHITECTURAL C ONTROL OF I CE -T EMPLATED C OLLAGEN S CAFFOLDS U SING A P REDICTIVE M ODEL

The architectural and physiomechanical properties of regenerative scaffolds have been shown to improve engineered tissue function at both a cellular and tissue level. The fabrication of regenerative three-dimensional scaffolds that precisely replicate the complex hierarchical structure of native tissue, however, remains a challenge. The aim of this work is therefore two-fold: i) demonstrate an innovative multidirectional freeze-casting system to afford precise architectural control of ice-templated collagen scaffolds; and ii) present a predictive simulation as an experimental design tool for bespoke scaffold architecture. We used embedded heat sources within the freeze-casting mold to manipulate the local thermal environment during solidification of ice-templated collagen scaffolds. The resultant scaffolds comprised complex and spatially varied lamellar orientations that correlated with the imposed thermal environment and could be readily controlled by varying the geometry and power of the heat sources. The complex macro-architecture did not interrupt the hierarchical features characteristic of ice-templated scaffolds, but pore orientation had a significant impact on the stiffness of resultant structures under compression. Furthermore, our finite element model (FEM) accurately predicted the thermal environment and illustrated the freezing front topography within the mold during solidification. The lamellar orientation of freeze-cast scaffolds was also predicted using thermal gradient vector direction immediately prior to phase change. In combination our FEM and bespoke freeze-casting system present an exciting opportunity for tailored architectural design of ice-templated regenerative scaffolds that mimic the complex hierarchical environment of the native extracellular matrix.


Introduction
The fabrication of regenerative scaffolds that replicate the structure and function of native tissue has proven a great challenge in recent years.An increasingly large body of literature has demonstrated the functional importance of organ-specific structure at multiple length scales [1][2][3][4][5][6][7].
The importance of extracellular architecture on tissue function is exemplified by many tissues and two examples are ventricular myocardium, which is composed of a transmural helical tissue that facilitates directed cardiac contraction; and osteons in bone, that are aligned longitudinally and comprise concentric helical layers that contribute to strength and toughness under directional load [8,9].In the context of regenerative tissue design, biomimetic structural and physiomechanical properties of regenerative scaffolds have been shown to improve cellular phenotyping [10], cellular function [11], native matrix deposition [12], and regenerated tissue function [13][14][15].The causal relationship between structure and function has heightened interest in development of controlled scaffold design and production techniques.
The programable nature of the technique affords advanced design and precise control of scaffold macrostructure, however, due to the extrusion process, scaffolds typically comprise of smooth solid struts that limit cellular adhesion, nutrient profusion, and long-range cellular communication [21,23,24].
Freeze-casting circumvents these limitations by leveraging the natural thermodynamics of the ice solidification process to define the multi-scale hierarchical architecture of resultant scaffolds [25][26][27][28][29].The low solubility of impurities in ice results in solute marginalization between ice grains during crystal formation [25,29].Thus, after sublimation, the resulting scaffold architecture is an exact negative of the ice crystal network [27,28].Manipulation of physical and chemical conditions during freezing has been shown to influence pore size [28,30,31], wall thickness [32], surface roughness [33], and interconnectivity [33,34].Additionally, by imposing a sufficiently high temperature gradient, directional crystal growth can be induced leading to structures with a preferred pore orientation [30,35].Some macrostructural diversity has also been demonstrated through ice templating, including isotropic [34,36,37], unidirectionally aligned [38][39][40], and radially aligned [40][41][42] pore morphologies.Further work also included ice nucleation from a tilted wedge [43,44] leading to bidirectional laminar structures, and from a radial geometry to radial and centrosymmetric growth [45].A strong relationship between thermal gradients and pore orientation has been established [25,27,46,47].To date, despite the well-defined thermodynamic process, no specific 3D control over complex lamellar orientation has been reported and macrostructural prediction and design remains a challenge.
The robust connection between scaffold architecture and the physics of ice growth has enabled numerous mathematical models and simulations of the freeze-casting process, including one dimensional particle diffusion models [48], phase field simulations [49,50], molecular dynamical simulations [51], and finite element models (FEM) [52].Many of these simulations are mathematically complex and computationally expensive making them difficult to use as an experimental design tool with the exception of FEM.Husmann et al. 2015 used a FEM to predict isotropic or anisotropic pore structure of resultant scaffolds based on the temperature distribution at the time of nucleation within various molds.This approach provided an accessible tool for experimental design.Due to the complexity of phase transition, however, no such tool has been created to predict spatially variable macrostructure in anisotropic ice-templated scaffolds [52].
In this work we present a bespoke multidirectional freeze-casting system, utilizing controlled heat sources integrated into casting molds to afford 3D thermal control over the solidification process and consequently precise control over lamellar orientation of resultant scaffolds.Directed control over ice templated collagen scaffold macro-architecture affords tissue specific design capabilities without compromising the micro-structural environment.Additionally, a finite element simulation was constructed to illustrate the influence of our system on the 3D thermal environment within the mold and to predict the lamellar orientation of resultant scaffolds.To test the efficacy of our approach collagen scaffolds were produced under three distinct thermal environments of increasing complexity.We measured time dependent temperature data and performed orientation analysis on resultant scaffolds to correlate the thermal environment during solidification with scaffold architecture.By varying the power output of embedded heat sources, we further assessed the capacity of our system to deliver adjustable and localized control of lamellar orientation.All of our findings were compared with simulation results to evaluate the utility of the FEM as an experimental design tool.Finally, we determine the impact of complex lamellar architecture on functional determinants of engineered tissue by characterizing the hierarchical structure and compressive mechanics of resultant scaffolds.

Freeze-casting system construction
Our multidirectional freeze-casting system was composed of a single controlled heat sink at the base of the casting mold and multiple heat sources embedded in the mold walls.The heat sink portion of the apparatus was constructed based on work presented by Deville et al. [26,53].However, in this work, only one thermally controlled heat sink was present at the base of the mold and hence the freezing front velocity was not dictated.Our freezing-system is shown in Figure 1A-B.In brief, the heat sink was controlled by a proportional-integral-derivative (PID) (Eurotherm 2400) controlled heater placed around a copper rod (2 cm Ø) cooled via liquid nitrogen with an ethanol buffer.
Our heat source mold system was composed of cylindrical polycarbonate molds composed of an external polycarbonate tube bound to a copper base and a removable nested internal polycarbonate tube (Figure 1C).Two resistance heaters (RS PT100 resistor) were embedded into the wall of the internal tube in the geometry depicted in Figure 1C.Heaters were individually controlled via a power supply (Tenma 72-10480 DC).

Freezing apparatus Efficacy
The thermal influence of the multidirectional freeze-casting system was tested by freezing water (9 ml) with the heat sink set to hold at −60 °C.The power input to the topmost heater (  ) was varied (  = 0, 0.06, 0.25, 0.56 and 1 W) for each iteration.

Finite element model construction
COMSOL Multiphysics 5.0 was used to create a finite element phase change simulation of the solidification process during freeze-casting.

Geometry and material properties
The simulation geometry and dimensions can be seen in Figure 1D and consisted of the threepart mold and a slurry cylinder.Two square regions (8x8 mm) on the external surface of the nested tube were partitioned to simulate the heat sources denoted   and   for the top and bottom respectively (Figure 1D).Thermal ties were imposed on all contacting surfaces via the form union function.The geometry was meshed with a physics defined 'finer' tetrahedral mesh, with a maximum and minimum element size of 1.5 and 0.1 mm respectively.Materials properties for polycarbonate and copper were imposed on their respective parts.The slurry cylinder was assigned the thermal and phase change properties of water as the thermal properties of collagen slurry have been shown to be comparable to water (Figure 1E-F) [52,54,55].
All material properties are summarized in Table 1.centered on the base of the mold.The heat transfer between the ambient air and geometry was simulated via convective heat flux, where   is the ambient temperature (24 °C),   is the temperature of the geometry surface, and ℎ  is the convective coefficient which is specific to the experimental setting [52,56].Values between 20 and 30 W m −2 K −1 incremented by 2.5 were simulated and the resultant steady state temperatures at the top of the slurry were compared to experimental data (section 3.1:   = 0).
Three thermal profiles were assessed; where heat source power input,   and   were varied.

Finite element model implementation
The 3D diffusion equation was used to calculate temperature, where T is temperature, t is time, Q is the sum of all heat sources in the system and ,   and k are material properties (Table 1).Due to the high viscosity of the slurry heat conduction via liquid convection was not considered.The phase transition was governed by equations: .
where  represents the number of phases and  represents the phase proportion and  is the phase change function: which has a rage of 0 to 1 for pure liquid and solid elements respectively.The phase transition was defined to occur over a 4 °C interval centered at 0 °C, based on previously reported differential scanning calorimetry data of 1 wt.% collagen suspension [37].Latent heat (333.5 kJ kg −1 ) was simulated via the apparent specific heat method and released during this interval [57].
Equations 2-7 were solved in conjunction with the boundary conditions presented in section 2.2.2 over a period of 3000 seconds at 1 second intervals.The final time-dependent temperature and temperature gradient solutions for all thermal conditions were exported for analysis.

Simulation analysis
All simulation analysis was done in Matlab® R2020a.

Thermal environment
Time dependent temperature data were averaged over 2 mm cubic regions within the slurry geometry.Regions were selected to correlate with thermocouple placement during experimental solidification.

Thermal gradient vector fields
The temperature gradient directions at the time of phase transition (time =  → = 0°C) were compared with experimental results.Gradient fields were visualized in cross section (−1 mm <  < 1 mm) as projections onto the xz-plane to replicate the planar slicing used for experimental analysis.Vector orientations were averaged over each of the nine regions of interest (ROIs) produced by dividing the x and z-axes into thirds as shown in Figure 1G.Standard deviation was used to indicate orientation variance within each ROI.

Directional freeze-casting
Suspensions of 1 wt% type I collagen (dermal bovine; Devro), were prepared in 0.05 M acetic acid and rested at 4 °C for 24 hours.Suspensions were homogenized at 22000 r.p.m. for 6 minutes and degassed via vacuum chamber (VirTis SP Scientific Wizard 2.0) ramped from 750 Torr to 2000 mTorr over 10 min.Directional freezing of the collagen slurry was executed with the freezing apparatus presented in section 2.1.1.Room temperature (~25 °C) suspension (9 ml) was pipetted into heat source molds and loaded onto a pre-cooled cold finger, held at −60 °C for all iterations.Heat source protocols varied for the three thermal conditions presented in section 2.2.1.
The single source condition was further modified by varying power output (  = 0.1, 0.2 and 0.65 W).Temperature data were collected (Omega RDXL6SD-USB data logger) for each iteration with k-type thermocouples at 4 second intervals.
Once freezing was completed the mold was immediately placed into a pre-cooled (−30 °C) freeze drier (VirTis SP Scientific Wizard 2.0) and dried at 0 °C under a vacuum of less than 100 mTorr for 20 hours.Lamellar orientation was assessed via Fourier power spectrum analysis with the Directionality plugin for ImageJ using the Fourier component method.The mean power for each lamellar angle was calculated across all slices within the ROI (Figure 1G).A Gaussian distribution was fit to the highest peak of the power spectrum.Directionality in the xz-plane for each ROI was identified as the peak of the Gaussian function and the standard deviation was calculated using Matlab® R2020a (Supp.Mat.S1-S2).The mean direction and standard deviation from corresponding ROIs of each scaffold reptation were compared.

Scanning electron microscopy (SEM)
Collagen scaffold sections were sampled from 2 cm above the base of the scaffold.These were sputter coated with gold for 2 min at a current of 20 mA and SEM micrographs were taken using a JEOL 820, with a tungsten source, operated at 10 kV.

Compression testing
Prior to compression testing, cross-linking was carried out using an EDC/NHS (1-ethyl-3-(3dimethylaminopropyl) carbo-diamide N-hydroxysuccinimide) chemistry (Fisher Scientific) using a molar ratio of 15:6:10 EDC:NHS:COOH (30% of the standard 5:2:1) [58,59].Reagents were dissolved in 95% ethanol and scaffolds were soaked for 2 hours at room temperature.Scaffolds were washed (5 x 5 minutes) with deionized water and freeze dried (VirTis SP Scientific Wizard 2.0) with a cooling rate of 0.2 °C min −1 to −20 °C.Drying occurred at 0 °C under a vacuum of less than 100 mTorr for 20 hours.Scaffolds were cut into cylinders (5 mm Ø; height: 5 mm) using a biopsy punch.All samples were punched at a height of 2 cm from the base of the scaffold.
Samples were composed of one of three architectures: 1. vertically aligned, 2. horizontally aligned, and 3. curved (35°6) alignment.1. and 2. were punched from scaffolds fabricated via the no source condition on the transverse and longitudinal planes respectively and 3. was punched from scaffolds fabricated via the single source condition, adjacent to the source on the transverse plane.
Hydrated (DI water) compression testing with a 5 kN Tinius Olsen machine and a 5 N load cell was executed.

Statistics
All the results are stated as mean ± standard deviation.Matlab® R2020a software was used for statistical analysis using a two-sample t-test where * indicates p<0.05.

Determination of heat transfer coefficient (𝒉 𝒄 )
The heat transfer coefficient (ℎ  ) is used in equation 1. to simulate convective heat transfer, it has been shown to lie between 5 and 30 W m −2 K −1 [52,56].A coefficient of ℎ  = 27.5 W m −2 K −1 most accurately described the experimental data (Table 2).

Thermal control achieved through multidirectional freeze-casting system
The localized impact of our heat source mold on the thermal environment during water solidification was evaluated.The temperatures across the height and width of the mold were measured (Figure 2).Thermocouple placement is shown in Figure 2A.The introduction of a heat source (  ) resulted in two independent temperature gradients (Figure 2B).The heat sink at the base of the mold (  = −60 °C) produced a vertical temperature difference due to the rapid temperature drop at the bottom of the slurry (Figure 2C).Powers above 0.1 W lead to a horizontal gradient and powers above 0.3 W to a change in vertical gradient (Figure 2B).A linear relationship between the horizontal temperature difference and   was identified ( = 25.1− 2.5; R 2 =0.995) (Figure 2B).Freezing times, especially adjacent to the source, were found to increase with power input such that powers above 0.6 W impaired solidification in this region (Figure 2D).

Simulation predictions and experimentally produced collagen scaffolds
We present our simulation and experimental results in tandem throughout the proceeding sections to highlight the agreement between our FEM simulation and the experimental freeze casting process.Initially, we compare the thermal environment during solidification under each source condition: no source (  =   = 0 ), single source (  = 0.56 ;   = 0 ), and Dual source (  =   = 0.56 ).We then assess the predictive capacity of the simulated temperature gradient at the time of freezing as an indicator for the final lamellar orientation of experimentally produced collagen scaffolds under each condition.Finally, we utilize the single source condition to explore heat source power output as an effective and predictable architectural control parameter.(Figure 3D).

Simulation predicts thermal environment during solidification
The dual source condition produced local temperature increases adjacent to both heat sources resulting in two horizontal gradients with antiparallel orientations in both the simulated and experimental results (Figure 3F & I).The temperature differences at   and   were simulated to be 15.1 and 15.0°C respectively (Figure 3E) and found experimentally to be 15.110.9 and 20.10.1°Crespectively (Figure 3F).Similarly to the no source condition, temperature elevations proximal to   were maintained, local temperature elevations near   , however, dissipated rapidly due to proximity to the heat sink (Figure 3I).Considering the continued presence of the vertical temperature gradient induced by the heat sink, the dual source condition resulted in three spatially independent thermal gradients during solidification Isothermal surface plots of simulated temperature results were used to describe the 3D thermal environment within the mold (Figure 3G-I).The no source simulation produced isothermal surfaces parallel to the base of the mold, indicating a uniformly perpendicular freezing front progression along the z-axis (Figure 3C).Introduction of localized heat sources distorted the freezing front, figures 3H-I the cooling fronts are shown to remain parallel to the base of the mold in regions distant from the heat sources, however the surfaces bend locally to become concave to the wall of the mold proximal to the sources.

Thermal gradient simulations at the time of phase change correlate with lamellar orientation in experimentally produced scaffolds
The gradient vector orientation across each mesh element at time  =  → was extracted from the simulation to create a 3D vector field prediction of final scaffold lamellar orientation.For comparison with scaffold architecture, the average gradient direction was calculated for each ROI as shown in Figure 4.The no source condition resulted in uniform orientations close to normal (90°) for both the simulation and experimental protocols (Figure 4A-C).Figure 4D-F shows that the introduction of   in the single source condition resulted in local deviations in both the gradient vector field and lamellar orientations.The mean orientation of the two ROIs (top and middle) adjacent to   was 121±13° and 125±6° for simulated gradients and scaffold lamellae respectively.Normal lamellar orientations, however, were both predicted and observed in regions distant from   (Figure 4D-F).
Simulation of the dual source condition resulted in a more complex vector pattern.The average orientations for the two ROIs adjacent to   and   (middle and bottom) were 120±15° and 76±10° respectively (Figure 4G).These predictions accurately described experimental findings where average lamellar orientations were 123±6° and 76±6° for the two ROIs proximal to   and   respectively (Figure 4H-I).
By comparing the resultant lamellar structure of scaffolds produced experimentally under each condition we found that the geometry of our multidirectional freeze-casting system significantly influenced lamellar orientation of resultant scaffolds (Figure 4J-M).The   source more effectively induced lamellar deviation from normal with 46 ± 2.7° and 40 ± 14.5° deviations observed in the top left ROI for the single source and dual source conditions respectively compared to 2±7 for the no source condition (p-values 3.9x10 -5 and 0.01) (Figure 4J).The middle left ROI was similarly affected by   (Figure 4K).In contrast the lamellar orientation shift that resulted from   had greater variability and reduced magnitude when compared to   .While the deviation from normal observed in the middle right ROI (−18.7 ± 11.0°) was significantly distinct from the same ROI produced under the single source condition ( = 0.04) (Figure 4L), no significant orientation shifts were observed in the bottom right ROI (Figure 4M).

Power output as a control parameter for lamellar orientation
Heat source power was assessed as a control parameter by varying the power output of   (  =0.1, 0.2, 0.56 and 0.65 W). Figure 5A shows simulated gradient vector fields for the two ROIs adjacent to   .These results accurately predicted the positive relationship between power output and lamellar deviation from normal as shown in Figure 5B.While both ROIs were influenced by   the top ROI was more sensitive to the heat source than the middle ROI in both the simulation results and experimental measurements Figure 5C.The power output also positively influenced the temperature difference across the width of the suspension and the time of complete phase transition (Table 3).Both parameters, however, were systematically underestimated in the simulation at all power outputs (Table 3).
There was a strong positive relationship between the horizontal temperature difference across the width (x-axis) of the mold at the time of freezing (∆( → )) and the resultant lamellar orientation for the single source geometry (Table 3).This pattern was not generalizable to sources at different heights within the mold.In the dual source condition, a temperature difference at phase change of 9.4±1.8°C at the height of   resulted in a lamellar deviation of 30±6° while a comparatively large temperature difference at freezing of 24.2±9.5 °C at the height of   resulted in a lamellar deviation of 14±6° (Figure 5D).

Lamellar orientation influences compressive stiffness
Compressive stiffness was measured for three lamellar compositions: parallel, perpendicular, and curved (average angle: 35±6°) relative to the loading force.The parallel and perpendicular samples were harvested from scaffolds produced under the no source condition.Samples with curved lamellar architecture were harvested from the region immediately adjacent to   of scaffolds produced under the single source condition.The modulus was found to be significantly increased for lamellar alignment parallel to the compressive load relative to the perpendicular and curved orientation (Figure 6C-D).Curved structures had intermediate stiffness that was comparable to structures with perpendicular orientation and was not statistically distinct (both pvalues 0.02; Figure 6D).

Efficacy of bespoke multidirectional freeze-casting system
In this work we developed a multidirectional freeze-casting system to afford 3D thermal control over the solidification process and consequently dictate the macrostructure of ice-templated collagen scaffolds.Previous studies have manipulated thermal sinks to alter gross pore geometry in ice-templated collagen scaffolds.These techniques, however, commonly rely on a single continuous heat sink and thus don't afford localized control over pore orientation [40][41][42]46].Our multidirectional freeze-casting system affords spatially specific macro-architectural control capabilities to more closely replicate native tissue without compromising the microstructural environment necessary for engineered tissue viability.By combing individually controlled heat sources and sinks we were able to effectively dictate the direction of ice propagation in a single step (Figure 2-Figure 3) [35].Resultant scaffolds had complex and spatially variable lamellar orientation dictated by directed processing parameters.Scaffolds consequently displayed variable intra-structural mechanical properties.Improved biomimetic architecture with hierarchical features will contribute to a wide range of disease modeling and clinical applications such as micro tissue, organoids, organ on a chip systems, surgical fixation devices and clinical regenerative implants [9,60].Additionally, due to the utilization of intrinsic solidification kinetics this method is scalable, however, because its efficacy relies on thermal gradients in space, the magnitude and geometry of heat sources and sinks must be scaled according to changes in mold dimensions.
The hierarchical architecture is due to the intrinsic crystallization process of ice.Anisotropic icetemplating has been shown to progress in four stages; nucleation, planar progression, columnar progression and lamellar/dendritic progression [30].Nucleation commonly occurs at a temperature below the equilibrium temperature and is followed by a rapid isotropic ice growth that results in a thin, dense isotropic region where solute is trapped in the planar solidification front [30].Ice front slowing due to synchronized latent heat release leads to a columnar progression and solute rejection from forming ice crystals.The thermal mass of the unfrozen suspension and the disequilibrium of solute concentration due to particle rejection further slows the ice front to a stable lamellar/dendritic progression.Ice crystals with preferential growth axes that sufficiently align with surrounding thermal gradients will grow at the expense of crystals with alternative orientations, and the true direction of ice propagation will fall between the direction of the external thermal gradient and the preferential growth axis of the [27,30] The intrinsic ice crystallization process, however, imposes limitations in reproducibility and scalability.Nucleation occurs below the equilibrium temperature and planar ice growth subsequently solidifies the supercooled suspension.Anisotropic scaffold size is therefore limited by the space needed to achieve sustained lamellar/dendritic ice front progression [30,61].Furthermore, reproducibility of ice templated scaffold characteristics such as surface roughness, microstructure and lamellar thickness has remained a challenge, due largely to the dendritic behavior of the lamellar growth phase in addition to stochastic nucleation and the isotropic nature of the planar phase [27].
The power output of embedded heat sources was as a control parameter for scaffold architecture.The positive relationship between lamellar deviation from the normal direction and heat source power demonstrates the adjustable control of local scaffold architecture (Figure 5C).
Power output directly influences the lateral temperature difference at the height of the source.This unidirectional parameter, however, is not a predictive indicator of lamellar orientation.The sensitivity of lamellar structure to lateral temperature difference was reduced at the bottom of the mold (  ) relative to the top (  ) (Figure 4Figure 5D).This can be attributed to the large vertical gradient component that dominates the thermal environment proximal to the heat sink.This result exemplifies the importance of a full 3D understanding of the thermal environment when manipulating the solidification process.For this reason the power output parameter is limited by the capacity for heat removal of the thermal sink; power output that exceeds this capacity will inhibit solidification and could damage the collagen protein (Figure 2) [46,62,63].Conversely, power output with a negligible magnitude relative to competing gradients will not influence crystal growth direction [46,47,64,65].Furthermore, the confounding relationship between power output and freezing kinetics necessitates consideration as there is a well-documented relationship between freezing kinetics and pore size [39,40,66].
We further characterized resultant scaffolds through SEM and compression testing in order to understand how our heat source protocol influenced functionality.One key strength of the icetemplating process is the inherent recapitulation of the hierarchical structure seen in the native ECM encompassing micro-and nanoscale features [21,25,42,67,68].These features have been shown to enhance cellular function by providing local environmental cues such as surface texture [21,[69][70][71] and improved nutrient profusion through microporosity [72][73][74].The surface topography and microporosity that comprised scaffolds produced under no source and single source conditions was consistent with conventionally freeze-cast structures [25,29,39] (Figure 6).Scaffold stiffness also an important functional modulator of engineered tissue [75][76][77][78].While cross-linking has been shown to increase the stiffness of collagen scaffolds it does not change the pore architecture, therefore enabling structural comparisons under load [79].Our results are consistent with composite mechanics where scaffold stiffness decreased rapidly as the angle of lamellar orientation deviated from parallel to the loading direction [80] (Figure 6).Taken together these results indicate that scaffolds that have locally variable lamellar alignment, such as those produced under the single source or dual source conditions, will have uniform microstructural features but distinct and locally variable intra-structural stiffness under load, a characteristic that has been documented in numerous native tissues [78].Scaffold microstructure and mechanical properties can be further tailored through additives or by varying collagen concentration which would allow finesse of scaffolds with bespoken properties [33,81,82].

Utility of predictive FEM as an experimental design tool
We presented our findings alongside a finite element simulation to explore the utility of a simple FEM as an experimental design tool.Previous FEM simulations of the ice-templating process have focused on determining the thermal environment at the time of nucleation and thus have not considered the phase change process [52].The present study aimed to describe the thermal environment during solidification of entirely anisotropic freeze-cast structures, necessitating the consideration of phase change within the simulation as the sequential release of latent heat during ice front progression contributes to the local thermal environment during solidification [49,66,83].
When compared to experimental data our simulation adequately described the time dependent temperature profiles of each condition, although the magnitude of the sustained lateral temperature differences, the time of phase transition ( → ) and the concurrent temperature difference (∆( → )) were systematically underestimated by the simulation (Figure 3;Table 3).
The simulated gradient vector field at   →  consistently predicted the lamellar orientation of experimentally produced collagen scaffolds for all ROIs (Figure 4).The simulation maintained predictive capacity when power output parameters were varied and accurately predicted the relationship between lamellar orientation and heat source power (Figure 5).This correlation indicates that the FEM accurately simulates the solidification process and substantiates the well documented correlation between temperature gradients and lamellar orientation [27,30,46].
Particularly, this agreement demonstrates the utility of our FEM as an experimental design tool in the fabrication of complex multidirectional freeze-cast structures.The simulation serves as a computational interface in which the influence of boundary conditions and geometries on the dynamic thermal environment within the mold during solidification can be explored and manipulated in advance to produce bespoke scaffold structures with tailored macro-architecture.Of note, it has been shown in [84,85] through very careful experiments, crystal orientation plays an important role in directional ice growth as ice has a strong preferred growth direction normal to its c-axis leading to tilted growth.In our experiments, we do not observe large deviations between pore alignments as seen in experiments compared to simulations (Figure 4).However, the increased orientation variance observed in the experimental scaffolds compared to the simulation could result from this competition between preferred directional growth and thermal gradient.
The simulation results, while largely in accordance with experimental findings, were subject to some limitations.A future experiment could look at the crystal orientation to explore this competition in more detail.Within the simulation, the phase change enthalpies were balanced using the apparent specific heat method, utilizing differential scanning calorimetry data of 1 wt.% collagen suspension reported by Pawelec et al. [37].The apparent specific heat method is computationally straightforward but has been shown to underestimate latent heat release and may have contributed to the underestimated time of phase transition produced by the simulation (Table 3) [57,86].The absence of nucleation stochastics in the model and the presence of solute in the experiment may also contribute to discrepancies in freezing kinetics between the experiment and simulation [18,28].Additionally, spatial discretization of data for analysis could have introduced confounding error as measurements were found to be sensitive to sampling location.
Having demonstrated the predictive power of the FEM, it is possible to use simulations to design the optimal experimental setup and achieve specific marco-architectures of resultant icetemplated scaffolds.The exciting future of this advance is to build a model which can take an ideal scaffold structure as an input and produce the experimental parameters (heater source powers, locations, temporal intervals) that will produce the desired architecture.

Conclusion
We have demonstrated the fabrication of hierarchically structured ice-templated collagen scaffolds with exceptional control of lamellar architecture.The heat sources embedded into freezecasting mold walls established and maintained a complex thermal environment during solidification.The lamellar orientation of resultant freeze-cast structures was dictated by the thermal environment during solidification and lamellar orientation was controlled locally by varying the source geometry and power output.Variations in macro-architecture influenced scaffold stiffness as curving lamellar structures significantly reduced the Youngs modulus under compression.Furthermore, we developed a predictive FEM to serve as an experimental design tool.Our simulation described the time dependent temperature profile of the cooling slurry and further illustrated the 3D thermal environment and freezing front topography within the mold during solidification.The simulated gradient vector direction at  → accurately predicted the lamellar orientation of freeze-cast collagen scaffold produced under analogous experimental conditions.By combining the predictive qualities of our FEM and the efficacy of our multidirectional freeze-casting system it is possible to tailor experimental design to create bespoke scaffolds with complex spatially varied lamellar architecture that recapitulates the native extracellular matrix and optimize regenerative tissue function and utility.

Figure 3
Figure 3 compares simulation and experimental time dependent temperature results from the

3. 4
Scaffold characterization 3.4.1 SEM micrographs show hierarchical scaffold structure SEM micrographs of scaffold samples from the no source condition were compared to samples from the region adjacent to   produced under the single source condition.Samples were examined for hierarchical pore structure and microscale surface features to determine if our heat source protocol disrupted scaffold microstructure.Both conditions resulted in microscale pores and surface features on the strut walls indicated by yellow and red arrows respectively in Figure 6A-B.

Table 2
Heat transfer coefficient estimation using the steady state temperature (SS temp) at the top of the slurry after solidification with the no source condition

Table 3
Influence of varied   on solidification parameters for the simulation and experimentally.
Solidification parameters include: deviation of the simulated thermal gradient angle and the lamellar orientation angle from the normal direction relative to the x,y-plane, time of solidification, and the temperature difference across the width of the mold