Towards a Comprehensive Framework for Made-to-Measure Alginate Scaffolds for Tissue Engineering Using Numerical Simulation

1. Introduction

1.1. Alginate Hydrogels

Alginate polysaccharides derived from brown algae are often used as excipients or bulk materials in biotechnological and biomedical products. In biomedical applications, these anionic polymers are applied as immobilization matrix of enzymes or cells, carriers for controlled release of pharmaceutical active compounds, wound dressings and tooth impression materials in dentistry.

In medical biotechnology, alginate hydrogels and their derivatives are used as scaffolds, microcarriers and immobilization matrices generating defined cell niches for physiologic cell-based models [1], e.g. for disease modeling [2]] or drug discovery [3]. Furthermore, alginate beads can be modified with extracellular matrix proteins and used in suspension bioreactors for the expansion of anchorage-dependent stem cells [4].

The versatile industrial use of alginates is attributed to their unique physico-chemical properties: they are non-toxic, biocompatible and gel in very mild conditions (ambient temperature, physiological pH, low concentrated and non-toxic cross-linking compounds).

Chemically, alginates are composed of unbranched polymers of the monomers β-D-mannuronate (M) and α-L-guluronate (G) which are linked via glycosidic bonds. The variation of the G/M ratio and the composition of homogeneous MMM- and GGG-polymer blocks or heterogeneous GMG/MGM polymer blocks are characteristic of different algae species and are influenced, e.g., by their local habitats or seasons. During cross-linking, mainly G-blocks interact with multivalent cations such as Ca²⁺, Ba²⁺, and Fe²⁺, forming the so-called egg-box structure. This interaction results in hydrophilic hydrogels with a water content of more than 98% [5]. The mechanical properties (e.g. stiffness, general stability) of gelled hydrogels depend on a variety of parameters such as the type of the cross-linking ion, its concentration, or overall cross-linking time. In the past decades, ultra-high viscosity (UHV) alginates [6,7] have been extensively studied and established as promising biomaterials in regenerative medicine and tissue engineering [8]. In contrast to other commercially available alginates, UHV-alginates differ mainly by having tenfold longer polymers while possessing the same chemistry and cross-linking mechanisms [8]. Furthermore, some special UHV-alginates are characterized by sterility, low endotoxin levels and overall high purity [9].

These special characteristics of the biomaterial result in low immunogenicity, high cell compatibility, high mechanical stability, higher chemical resistance, and higher elastic modulus, making this biomaterial suitable for applications in medical biotechnology.

1.2. Simulation/Model of Alginate Hydrogel Formation

Throughout the past few decades, simulations have provided a valuable and insightful way to gain insight into the formation of alginate hydrogels. For simplification, these mathematical models tend to focus on a few, or sometimes a single, aspect of this process. Stewart et al. [10] created a molecular dynamics (MD) simulation to investigate the complexation of poly-M and poly-GM oligomers with sodium and calcium ions, resulting in an alginate gel. Previous studies focused primarily on poly-G sequences [11,12,13]. They observed that all conducted simulations resulted in complexation with Na⁺ and Ca²⁺ ions, regardless of the used oligomer sequence. The formation of the „egg-box” structure appeared to follow the order of GG > MM > GM in terms of reaction preference. These findings suggest that poly-M regions may play a role in secondary cross-linking interactions.

Nair et al. [14] determined the mechanical properties of alginate-based 3D scaffolds using an in silico model, created from experimental data. They found that the Ogden hyperelastic model was the most suitable to describe the nonlinear behavior of alginates, compared to other hyperelastic material models. This model was then modified, so that the material constants became a function of alginate concentration and viscosity, which allowed them to analyze the effects of physical and geometric parameters such as porosity on the stress-strain curve. All their tests were validated against experimental data and showed a good level of agreement.

A micro-mechanical constitutive model was developed by Külcü et al. [15] to describe the nonlinear inelastic properties of alginate double network hydrogels by characterizing the elongation of the polysaccharide polymer chain after it has been fully stretched. They employed the Arruda-Boyce eight-chain model [16], the network decomposition approach [17,18] and the network alteration concept [19] to create a new energy function, describing the enthalpic and entropic behavior of a single polymer chain. This approach yielded good agreement with experimental results.

Mikkelsen et al. [20] proposed a mathematical model that describes calcium-induced 3D alginate gel formation using a set of coupled chemical reaction-diffusion equations. They used measurements from 1D gels to provide the necessary parameter values. Thu et al. [21] took experimentally determined gradients, gathered through MR-microimaging and SRIXE, from 3D gels to reevaluate these previous parameter estimates, e.g. the relative alginate diffusion coefficient and the reaction rate constant. The updated simulation model was then capable of producing results in good agreement with both 3D and 1D experimental data, obtained under various experimental conditions. Data comparisons supporting these claims were not included. Bjørnøy et al. presented an experimental framework for the gelation of low viscosity alginate with calcium concentrations higher than 50 mM. The main focus of their work is on the experimental toolbox to understand the gelation and mineralization of alginate hydrogels, only minor emphasis placed on the numerical modeling of this process [22].

1.3. Analysis of Hydrogel Cross-Linking Process

So far, numerous in vitro studies have been conducted to understand the correlation between type and concentration of the cross-linker and the resulting hydrogels. Since it is very challenging to differentiate between the transition of sol and gel states of hydrogels, researchers have found several approaches to quantify this process. Huynh et al. studied the gelation kinetics of polygalacturonate with four different divalent cations by using the changes of turbidity during sol-gel transition. With this method, they were able to determine, for examples, the diffusion coefficients of different cross-linker concentrations [23]. Despite the promising results of this method, it depends on the changing turbidity during gelation and is thus not universally applicable to hydrogels. Another approach for visual analysis of ion occurrence was presented by Cheng et al. In the context of an electro-induced gelation mechanism, they utilized a pH-dependent release of Ca²⁺ ions from CaCO₃ and monitored the gelation process with a Ca²⁺-dependent fluorescent dye. Although this method delivers data to follow the alginate cross-linking kinetics, it is primarily limited to Ca²⁺ induced gelation processes due to the specificity of available ion or metal indicators [24,25]. Photon Correlation Imaging was employed to investigate the cross-linking process in small volumes of alginates by Secchi et al. [26]. Hajikhani et al. proposed a chemo-mechanical model for alginate cross-linking based on their experimental data with commercial sodium alginate for bioprinting applications. They focused on diffusion and reaction kinetics, incorporating the resulting finite strains and shrinkage effects caused by syneresis. Additionally, they accounted for the crosslinking reaction rate dependence on internal gelling mechanics, establishing a two-way coupling between chemistry and mechanics [27,28]. Stößlein et al. proposed a method that allows the in-situ determination of alginate cross-linking by mechanical texture analysis. In their study, they compressed large alginate beads until rupture after different cross-linking times and measured the applied forces. They demonstrated a strong correlation between guluronate content and applied pressure of different types of alginates [29]. Cesewski et al. described in their study a cantilever-based method for real-time monitoring of alginate cross-linking process [30], while Potter et al. studied the cross-linking of alginate with Ca²⁺ process by magnetic resonance imaging [31]. Bjørnøy et al. analyzed the cross-linking kinetics of low-viscosity alginates in the presence of Ca²⁺. In their work, they imaged the cross-linking front by dark field and confocal laser scanning microscopy of regularly formed alginate discs or beads [32]. In their study, they used a modified numerical model originally presented by Thu et al. [21]. Furthermore, they studied only low- and medium- viscosity alginates and used cross-linker concentrations higher than 50 mM [32]. Hu et al. investigated the diffusion properties of various alginate formulations using different gelation baths. However, they performed experiments and tested the diffusion properties of the resulting beads with fluorophore-labeled bovine serum albumin but did not conduct any numerical simulation [33]. Li et al. [34] studied the gelation of alginate using different cross-linker concentrations and measured, for example, the swelling behavior, transparency, or FT-IR spectra. However, despite the new insights gained in alginate interaction with Ca²⁺, only low or medium-viscosity alginates were used, and numerical modeling has not been performed.

Besiri et al. used a rheometer approach with a modified rheometric setup to analyze the gelation of technical low viscous alginate using CaCl₂ [35].

An assessment of the polymer gradients using synchrotron radiation-induced x-ray emission, magnetic resonance microimaging, and mathematical modeling [21] was conducted by Thu et al. The study primarily focused on the inhomogeneity of alginate microbeads and not on the gelation kinetics. Furthermore, ultra-high-viscosity alginates were not considered in the study and gelation conducted with Ca²⁺-ions.

The numerical model presented in this study employs a reaction-diffusion approach to the gelation process. While others [22,28] have taken, modified and expanded on the framework developed by Mikkelsen and Elgsaeter [20], this work aims to explore a different numerical approach. Additionally, this study focuses on barium instead of calcium as the crosslinking ion in the process of egg-box formation.

1.4. Alginate Hydrogels in Tissue Engineering and Regenerative Medicine

In biotechnological applications, alginate hydrogels have been used for several decades primarily as an immobilization matrix for therapeutically relevant cells [36,37,38,39]. The discovery of induced pluripotency [40,41] with the possibility of deriving virtually any somatic cell type of the human body in vitro was a game-changer for regenerative medicine and tissue engineering. With the availability of these cell systems, there has been a growing demand for cell type-specific environments that replicate in vivo conditions more accurately compared to standard plastic dishes [42,43,44].

The rational design of the stem cell niche to provide tissue-specific environments combines soluble factors, adjacent cell types, specific proteins of the extracellular matrix, the topography of the adhesion site, and the mechanical properties of the bulk material. Cell-type adapted niches lead to in vivo-like behavior of the cells and, in consequence, to valuable data when these cells are used as models for drug screening, cytotoxicity and disease modelling. If these properties differ from the native environment, e.g. the cultivation on (stiff) 2D surfaces instead of the actual 3D environment [42,45], it could consequently lead to incorrect dose responses if these models are used in drug testing screenings, e.g., for personalized medicine, as shown for breast cancer drugs by Imamura et al. [46].

The mechanical properties of native tissue vary from very soft (e.g., brain) to ultra -stiff (e.g., bone) [47]. In alginates, gelation is the fundamental mechanism responsible for the mechanical properties. As described in section 1.1. multivalent cations interact with adjacent guluronic acid blocks and form the hydrogel. Depending on the chemical composition of the alginates used as well as the reaction conditions (time, cation concentration), softer or stiffer hydrogels can be produced. The alginate hydrogels can be used as a soft growth surface for stem cells [1,4] or as immobilization matrix for therapeutically relevant cells [8]. The mechanical stimuli, either from the surface or the surroundings, are linked through cell-matrix contacts to the interior of the cell and can trigger molecular pathways (mechanotransduction) [48,49]. In pluripotent stem cells, the so-called rigidity-dependent YAP/TAZ signaling pathway influences the cell during neural differentiation [50].

1.5. Rationale Predicting Alginate Hydrogel Gelation

In cell biology, protocols for generating soft 2D and 3D environments are available but are primarily driven by concentration and time. Furthermore, the final properties of the hydrogels are often characterized empirically, and changes in gelation conditions can only lead to a relative characterization of the hydrogel (“softer”, “stiffer”). Besides the chemical conditions of the reaction, the volume of the alginate to be gelled can be altered, e.g., if the cell culture substrate must be changed. Higher volumes with the same gelling conditions lead to inhomogeneous or only partailly gelled hydrogels. In general, the produced hydrogels must be reliable and their mechanical properties comparable, but an intensive quality control and characterization is often not feasible, require special hardware and expertise, and is in addition time consuming.

Gelation is also affected by encapsulated cells, whether single cells or multicellular aggregates. In the worst case, this influence can lead to critical structures or inhomogeneous gelation, which potentially triggering biological responses. A possible example of this is immuno-transplantation [9]. Overall, a predictive tool for the gelation of hydrogels is in high demand and would also be useful for cell biologists producing tissue-specific hydrogel surfaces or for quality control in regenerative medicine. In this study, the fundamental gelation kinetics will be modeled based on experimental data and validated in initial real-world scenarios.

The gelation history also significantly impacts the resulting properties of a cross-linked alginate construct [51]. For example, the diffusion properties of cross-linkers change noticeably after the alginate solution undergoes the cross-linking process [52]. This also affects the availability of open alginate cross-linking sites [53]. Both of these factors may result in irregular cross-linking within a gelled alginate structure, as well as mobility, concentration, and mechanical gradients, which conversely impact local cross-linking rates [51].

Overall, a numerical model that can predict the gelation behavior and the mechanical properties of alginate hydrogels is still lacking. Furthermore, a tool that can assist scientists to gel alginates for distinct biological cell models is currently unavailable. In this work, the fundamental numerical simulation for a software tool will be presented. This tool will be further enhanced in the future to provide biologists and material scientists with a virtual supervisor to produce highly defined alginate hydrogels for their applications.

2. Materials and Methods

2.1. Disc-Shaped Alginate Gelation Kinetics Study

Alginate hydrogels from brown algae Lessonia nigrescens and Lessonia trabeculata were purchased from alginatec GmbH (Riedenheim, Germany). In this study, a 1:1 mixture (vol%/vol%) of both alginate species were used. The alginates were solved in saline (0.9%, B. Braun Melsungen AG, Melsungen, Germany). The dynamic viscosity of the mixture was characterized using the rheometer MCR 92 (Anton Paar GmbH, Ostfildern, Germany) and was 3306 (± 78 mPa∙s). The different cross-linking solutions were prepared with 115 mM NaCl and 10, 20, and 40 mM BaCl₂ (Merck KGaA, Darmstadt, Germany). The osmolarity of all cross-linking solutions were adjusted to 270 mOsm/l. Alginate sol were pipetted on Nunc EasYDish dishes (ThermoFisher Scientific GmbH, Dreieich, Germany) and covered with silicone spacer (Grace Bio Labs, Merck KGaA, Darmstadt, Germany, CoverWell incubation chamber, 20 μl volume, 13 mm diameter and 200μm height). To facilitate the injection of cross-linking solutions, two thin channels were punched into the silicone wall. Dishes and silicone spacer were cleaned with desalted water in an ultrasonic bath and air dried before usage. One droplet of alginate (4 μl) has been carefully pipetted in the middle of the petri dish and covered with the silicone spacer. Enclosed air bubbles occurred rarely and did not influence the gelation process. Image acquisition has been started and 20 μl of the cross-linker solution has been pipetted through the inlet to fill the whole chamber and start the gelation process. Time lapse movies were acquired using a standard phase contrast microscope (Eclipse Ts2, Nikon Instruments, Duesseldorf, Germany) with 4x magnification. Images were captured automatically with a mounted camera (DMK 33UX174 and software IC capture 2.4, both The Imaging Source Europe GmbH, Bremen, Germany) in the following intervals: 0 to 5 min: each 10 s, 5 min to 15 min: each 20 s and for the rest: each 40 s. Captured images were analyzed using the IC Measure software (The Imaging Source Europe GmbH, Bremen, Germany).

Figure 1. Setup and principle of observing the gelation process. (a) A thin disc-like volume of alginate is poured into a dish and covered by a thin silicone spacer for gelation with different BaCl₂-solutions. This process can be observed using phase contrast microscopy and a concentric decrease in the travelling liquid/gelled interface can be tracked and used for quantification of the gelation process. (b) Schematic drawing in two different time points of alginate gelation. The disc-like volume of alginate is surrounded by the BaCl₂ cross-linker solutions and, consequently, barium (and chloride) ions diffuse into the alginate sol triggering the gelation that can be tracked by the traveling liquid/gelled interface over time. The diameters of the circular interfaces decrease over time and disappear after complete gelation of the alginate discs. (b) generated with BioRender.com.

Figure 1. Setup and principle of observing the gelation process. (a) A thin disc-like volume of alginate is poured into a dish and covered by a thin silicone spacer for gelation with different BaCl₂-solutions. This process can be observed using phase contrast microscopy and a concentric decrease in the travelling liquid/gelled interface can be tracked and used for quantification of the gelation process. (b) Schematic drawing in two different time points of alginate gelation. The disc-like volume of alginate is surrounded by the BaCl₂ cross-linker solutions and, consequently, barium (and chloride) ions diffuse into the alginate sol triggering the gelation that can be tracked by the traveling liquid/gelled interface over time. The diameters of the circular interfaces decrease over time and disappear after complete gelation of the alginate discs. (b) generated with BioRender.com.

2.1. Sphere-Shaped Alginate Gelation

Alginate hydrogels from the brown algae Lessonia nigrescens and Lessonia trabeculata were purchased from alginatec GmbH (Riedenheim, Germany). For the disc-shaped alginates, a 1:1 mixture (vol%/vol%) of both alginate species were used. The alginates were dissolved in saline (0.9%, B. Braun Melsungen AG, Melsungen, Germany). The dynamic viscosity of the mixture was characterized using a rheometer (MCR 102, Anton Paar GmbH, Ostfildern, Germany) and was 2041 (± 12) mPa∙s.

To determine the crosslinking kinetics of alginate beads, different crosslinking solutions of barium chloride (Merck KGaA, Darmstadt, Germany) were prepared at concentrations of 10 mM, 20 mM, and 40 mM. The osmolarity ranged from 291 to 301 mOsmol/l. 100 µL of each crosslinking solution was transferred into a 96 U-bottom well plate (Greiner Bio-One GmbH, Frickenhausen, Germany). The addition of alginate to the crosslinking solution was performed by dropping 5 µL with a pipette. Video sequences of the crosslinking process were captured using a 4x phase contrast objective. For observation of shrinkage, additional images were acquired after 180 min, while the samples remained in their respective barium chloride solutions. A Nikon Eclipse microscope TS2 coupled with a DS-Ri2 camera (Nikon Instruments, Duesseldorf, Germany) was used for this analysis. The samples were kept at room temperature. Video sequences and images were processed using NIS-Elements Advanced Research (NIS AR 5.41.00, Nikon Instruments, Duesseldorf, Germany).

The cross-linking progress was analyzed every five seconds, starting from 15 to 20 seconds after dropping the alginate into the cross-linking bath. The cross-linking progress was observed as long as the cross-linking front was visible, for a maximum of 290 seconds. The diameters of the cross-linking fronts were plotted for all five replicates over time, and the speed of the cross-linking front was determined by linear fit.

2.2. Modelling

2.2.1. Basic Assumptions

In the following investigation it is assumed that the alginate hydrogel consists of a polymer ${\phi }^{poly}$, an aqueous solvent ${\phi }^w$ and barium ions ${\phi }^{ba}$. During the cross-linking the barium ions form bonds between the molecules of the polymer resulting in the cross-linked polymer ${\phi }^{cross}$ (see also Figure 2a). The model is developed on the basis of a multi-phase continuum description, i.e. a mixture theory is chosen as starting point. Therefore, not the individual molecules or ions are investigated but certain amounts of them occupying a representative volume element.

The amount of mass of the individual molecules is represented by the partial densities ${\rho }^{\alpha }$ where $\alpha \in$ {poly, w, ba, cross} indicates the above-mentioned constituents.

Figure 2b shows a representative volume of size (RVE) $dv$ which contains $N^{poly}$ mols of polymer molecules each of molar mass $M^{poly}$, $N^{ba}$ free barium ions of molar mass $M^{ba}$ and $N^w$ water molecules of molar weight $M^w$. Only one cross-linked polymer molecule is present in the RVE. Its molar mass is increasing during the cross-linking according to the number of polymers being added to the cross-linked one and the required number of barium ions. In consequence the number of polymers and of barium ions is decreasing during cross-linking. The water content is not influenced by the cross-linking because it is acting as a solvent only.

The partial density of the constituents inside the RVE is defined as the ratio of its mass $d{m}^{\alpha }$ to the occupied volume $dv$. Using the number of molecules per volume $N^{\alpha }$ and the molar weight $M^{\alpha }$ the partial density of the constituent ${\phi }^{\alpha }$ in the RVE is expressed as

$${\rho }^{\alpha }={\displaystyle \frac{d{m}^{\alpha }}{dv}}={\displaystyle \frac{N^{\alpha }{M}^{\alpha }}{dv}}$$

(1)

While the total mass in the RVE is equal to the sum of the partial masses the mixture density is obtained as

$$\rho =\sum _{\alpha }{\rho }^{\alpha }.$$

(2)

2.2.2. Mixture Theory

In the following a continuum mixture theory is applied as the basis of the model. The concept assumes that the RVE introduced in section 2.3 is small compared to the macroscopic dimensions and, therefore, it is treated as a spatial point in the continuum description located at a position x. This leads to the concept of superimposed continua because the RVE is simultaneously filled by all constituents, each of them described by its partial mass density.

In addition to the mass exchange between the constituents ${\phi }^{poly}$, ${\phi }^{ba}$ and ${\phi }^{cross}$ the mass in the RVE can also change due to the flux of molecules/ions of the individual constituents over the boundary of the RVE. Balancing the mass flux ${\rho }^{\alpha }{v}^{\alpha }$ over the boundary of the RVE with the change of mass $\partial {\rho }^{\alpha }/\partial t$ inside the volume and the mass exchange ${\widehat{\rho }}^{\alpha }$ leads to the local formulation of the mass balance of the constituent ${\phi }^{\alpha }$

$${\displaystyle \frac{\partial {\rho }^{\alpha }}{\partial t}}+div({\rho }^{\alpha }{v}^{\alpha })={\widehat{\rho }}^{\alpha }.$$

(3)

Here, $v^{\alpha }$ is the macroscopic velocity of the constituent ${\phi }^{\alpha }$ and the operator div is the divergence with respect to the spatial position x. The situation is captured in Figure 2b, where a magnification of the point at location x shows the RVE.

Following Truesdell’s metaphysical principles of mixture theory the balance of total mass is obtained by summing up the mass balances of the individual constituents. Furthermore, the format of this balance must be identical to the balance equation of a single-phase continuum. Summation of (3) yields

$${\displaystyle \frac{\partial \sum _{\alpha }{\rho }^{\alpha }}{\partial t}}+div\sum _{\alpha }{\rho }^{\alpha }{v}^{\alpha }=\sum _{\alpha }{\widehat{\rho }}^{\alpha }.$$

(4)

Figure 2. (a) Representative volume filled with the free polymer, barium ions and cross-linked polymer and water (not shown). (b) Macroscopic domain and RVE as magnification of a spatial point. The mass of constituent φ^α inside the RVE changes due to the flux over the boundary and the mass exchange. Created with BioRender.com.

Figure 2. (a) Representative volume filled with the free polymer, barium ions and cross-linked polymer and water (not shown). (b) Macroscopic domain and RVE as magnification of a spatial point. The mass of constituent φ^α inside the RVE changes due to the flux over the boundary and the mass exchange. Created with BioRender.com.

This result is identical with the mass balance of a single-phase continuum if (2) is substituted into the first term and if the barycentric velocity or mixture velocity is defined according to

$$\rho v=\sum _{\alpha }{\rho }^{\alpha }{v}^{\alpha }.$$

(5)

Furthermore, the total mass is conserved. This requires that the sum of the mass exchange terms is zero

$$\sum _{\alpha }{\widehat{\rho }}^{\alpha }=0.$$

(6)

In the same manner further balance equations can be derived, e.g. the balances of momentum and of energy, cf. e.g. Truesdell [54], Bowen [55], etc. As a simplification it is assumed that during the cross-linking process the barycentric velocity is zero, i.e. there is no average macroscopic motion. But the individual constituents may move. To describe this relative motion of ${\phi }^{\alpha }$ with respect to the center of mass of the mixture, which stays at rest, the diffusion velocities are introduced by

$$d^{\alpha }=v^{\alpha }-v.$$

(7)

Using the definition (2) of the mixture mass density and the definition (5) of the barycentric velocity it can be proven that

$$\sum _{\alpha }{\rho }^{\alpha }{d}^{\alpha }=0.$$

(8)

With the definition of the barycentric velocity (5) and the diffusion velocity (7) the mass balances can be re-written as

$${\displaystyle \frac{\partial {\rho }^{\alpha }}{\partial t}}+div\left({\rho }^{\alpha }v+{\rho }^{\alpha }{d}^{\alpha }\right)={\widehat{\rho }}^{\alpha }.$$

(9)

In further investigations it is assumed that the diffusion velocities and their changes are small. Therefore, the associated change of momentum is small. This avoids a detailed study of the balances of momentum for the individual constituents and the diffusive mass fluxes

$$j^{\alpha }:={\rho }^{\alpha }{d}^{\alpha }$$

(10)

can be prescribed by constitutive equations. In the present case diffusion is the main driving force. Therefore, the flux terms are governed by Fick’s second law. For simplicity the barycentric motion is not further investigated during the cross-linking, i.e. no external loads are applied to the mixture,

$$v=0.$$

(11)

Based on the above-mentioned assumptions the cross-linking model consists of the four balance of mass equations

$${\displaystyle \frac{\partial {\rho }^{poly}}{\partial t}}=-div{j}^{poly}+{\widehat{\rho }}^{poly},$$

(12)

$${\displaystyle \frac{\partial {\rho }^w}{\partial t}}=-div{j}^w+{\widehat{\rho }}^w,$$

(13)

$${\displaystyle \frac{\partial {\rho }^{ba}}{\partial t}}=-div{j}^{ba}+{\widehat{\rho }}^{ba},$$

(14)

$${\displaystyle \frac{\partial {\rho }^{cross}}{\partial t}}=-div{j}^{cross}+{\widehat{\rho }}^{cross}$$

(15)

and the constraints

$$\sum _{\alpha }{\widehat{\rho }}^{\alpha }=0,$$

(16)

$$\sum _{\alpha }{\rho }^{\alpha }{d}^{\alpha }=0.$$

(17)

2.2.3. Constitutive Equations

For the following investigations the aqueous suspension is neglected, and the polymer, cross-linked polymer and barium constituents are investigated only. During the cross-linking the mass density of the cross-linked polymer is increasing while the densities of the barium ions and the polymers are decreasing. Considering the constraint (6), the mass production for the cross-linked polymer is given by

$${\widehat{\rho }}^{cross}=-{\widehat{\rho }}^{poly}-{\widehat{\rho }}^{ba}.$$

(18)

The probability for the formation of a new node in the cross-linked polymer network is increasing with the amount of available barium ions and with the amount of free polymers. These effects are reflected by the constitutive assumptions

$${\widehat{\rho }}^{poly}=-C^{poly}{\rho }^{poly}{\rho }^{ba},$$

(19)

$${\widehat{\rho }}^{ba}=-C^{poly}{\rho }^{poly}{\rho }^{ba}-C^{cross}\left(1-\kappa \right){\rho }^{cross}{\rho }^{ba}.$$

(20)

While one polymer can link several ions to the network, a second term is introduced in (20) with the following interpretation: While $C^{poly}$ describes the initial connection of one polymer to the cross-linked polymer, which requires one barium ion, a further cross-linking can take place. After the polymer became part of the cross-linked polymer, the remaining free positions can absorb further ions until all possible positions are occupied by barium ions and the maximum amount of cross-links is formed in the cross-linked polymer. This effect is described by the additional internal variable $\kappa$ defined as the ratio of occupied bindings to possible bindings. The term $\left(1-\kappa \right){\rho }^{cross}$ is a measure of the free binding positions. The evolution of κ is driven by a reaction equation

$${\displaystyle \frac{\partial \kappa }{\partial t}}=C^{\kappa }\left(1-\kappa \right){\rho }^{cross}{\rho }^{ba}.$$

(21)

During the formation of the polymer network, the diffusion motion of the cross-linked polymer is neglected, i.e.

$$d^{cross}=0.$$

(22)

This assumption implies that a swelling of the cross-linked polymer is currently not included in the model for simplicity. In contrast the polymer constituent and the barium ions move with respect to the stationary cross-linked polymer. When the motion is mainly driven by diffusion an appropriate choice for the mass fluxes is Fick’s law

$${\rho }^{poly}{d}^{poly}=-D^{poly}grad{\rho }^{poly},$$

(23)

$${\rho }^{ba}{d}^{ba}=-D^{ba}grad{\rho }^{ba}.$$

(24)

The size of ions is much smaller than the size of the polymers, therefore, the diffusion constant $D^{ba}$ is several orders of magnitude bigger than $D^{poly}$

$$D^{ba}\gg D^{poly}.$$

(25)

During the cross-linking the mesh size of the cross-linked polymer is decreasing. This drastically hinders the diffusion motion of the polymers but also the diffusion motion of the barium ions. For both diffusing constituents an exponential function is chosen relating the current value $D^{\alpha }$ of the diffusion constant to its initial value $D_0^{\alpha }$

$$D^{\alpha }=D_0^{\alpha }\exp \left(-K^{\alpha }{\rho }^{cross}\right)for\alpha \in \left\{poly,ba\right\}.$$

(26)

Due to the different sizes of the ions and the polymers the constants $K^{\alpha }$ are different.

2.3. Implementation

A set of kinetic differential equations based on the previous chapter is used to describe the local accumulation and depletion of mobile alginate, cross-linked alginate and cross-linking cations [54]. These kinetic differential equations are given as a set of diffusion-reaction equations (27-30).

$$\begin{array}{c}{\displaystyle \frac{\partial }{\partial t}}{\rho }^{poly}\left(t,\mathit{x}\right)=div\left(D_0^{poly}\exp \left(-K^{poly}{\rho }^{cross}\left(t,\mathit{x}\right)\right)div{\rho }^{poly}\left(t,\mathit{x}\right)\right)\\ -C^{poly}{\rho }^{poly}\left(t,\mathit{x}\right){\rho }^{ba}\left(t,\mathit{x}\right)\end{array}$$

(27)

$${\displaystyle \frac{\partial }{\partial t}}{\rho }^{ba}\left(t,\mathit{x}\right)=div\left(D_0^{ba}exp\left(-K^{ba}{\rho }^{cross}\left(t,\mathit{x}\right)\right)div{\rho }^{ba}\left(t,\mathit{x}\right)\right)$$

$$-C^{ba}{\rho }^{poly}\left(t,\mathit{x}\right){\rho }^{ba}\left(t,\mathit{x}\right)$$

$$\begin{array}{c}-C_{\kappa }^{ba}{\rho }^{ba}\left(t,\mathit{x}\right)\left(1-\kappa \left(t,\mathit{x}\right)\right){\rho }^{cross}\left(t,\mathit{x}\right)\end{array}$$

(28)

$${\displaystyle \frac{\partial }{\partial t}}{\rho }^{cross}\left(t,\mathit{x}\right)=\left(C^{poly}+C^{ba}\right){\rho }^{poly}\left(t,\mathit{x}\right){\rho }^{ba}\left(t,\mathit{x}\right)$$

$$\begin{array}{c}+C_{\kappa }^{ba}{\rho }^{ba}\left(t,\mathit{x}\right)\left(1-\kappa \left(t,\mathit{x}\right)\right){\rho }^{cross}\left(t,\mathit{x}\right)\end{array}$$

(29)

$${\displaystyle \frac{\partial }{\partial t}}\kappa \left(t,\mathit{x}\right)=C^{\kappa }{\rho }^{ba}\left(t,\mathit{x}\right)\left(1-\kappa \left(t,\mathit{x}\right)\right){\rho }^{cross}\left(t,\mathit{x}\right)$$

(30)

This numerical model was implemented using the FEM software COMSOL Multiphysics® v6.2. (COMSOL AB, Stockholm, Sweden) The simulations are run on a workstation containing a 16-core 11th Gen Intel® Core™ i9-11900K CPU and 64 GB of RAM.

The model parameters, which are given in Table 1, were optimized to achieve the best fit between the simulation results and the experimental data.

3. Results

3.1. Alginate Gelation

The gelation process from the alginate sol to the alginate hydrogel was analyzed by tracking a thin traveling front line from the outer edge to the center of the alginate disc. In a non-invasive analytical approach, this process was quantified using phase contrast microscopy and subsequent manual image analysis. Representative time-lapse sequences of cross-linking experiments are illustrated in Figure 3a (10 mM BaCl₂), Figure 3b (20 mM BaCl₂), and Figure 3c (40 mM BaCl₂). Due to the low contrast of the gelation interface, images in these figures were enhanced for better visibility by FFT Bandpass Filter of ImageJ (Fiji, parameters “filter large structures down to” 50 pixels, filter small structure up to 3 pixels, tolerance of direction 5%, with activated autoscaling and saturate image option).

The data reveal that the general course of the traveling liquid/gelled interface is comparable for all studied cross-linker concentrations (10 mM, 20 mM, and 40 mM, Figure 4a). The kinetics of all studied concentrations can be classified into two phases: the first phase at the beginning with a fast-traveling visible front line followed by a phase with a slower but also linear course. In addition, we observed a more decreasing slope/velocity in Figure 4a with increased cross-linker concentration, resulting in a faster gelation time of the alginate drop (Figure 4b). The velocity of gelation in this work is defined as reduction of ungelled core and negative.

The investigation of the cross-linking kinetics of three different BaCl₂ concentrations with comparable osmolarities was performed using UHV-alginate microcarriers by dropping 5 µL into the cross-linking bath. The cross-linking front was observed microscopically, and the advance of the front was measured by analyzing the diameter of the front over time. Since the sizes of the spheres differ from 1210 to 1705 µm, the initial sizes of the cross-linking fronts differ as well and require accordingly more or less time to reach the core of the spheres. The cross-linking fronts in the replicates with 10 mM BaCl₂ need 145 to 290 seconds to reach the core and are strongly dependent on the sphere’s diameter.

Figure 5b shows the diameter of the cross-linking front over time for a 20 mM BaCl2 cross-linking concentration. In all five replicates, the decrease in size over time was observed, with a steeper decrease than in Figure 5a. Since the sizes of the spheres differ from 1124 to 1885 µm in size the initial size of the cross-linking front differs as well and needs accordingly more or less time to reach the core of the spheres. The cross-linking fronts in replicates with 20 mM BaCl₂ require 100 to 200 seconds to reach the core and are strongly dependent on the sphere’s diameter.

Figure 5c shows the diameter of the cross-linking front over time for 40 mM BaCl₂ cross-linking concentration. In all five replicates, the decrease in size over time was observed, with the steepest decrease compared to Figure 5a and Figure 5b. The sizes of the spheres have comparable diameters (1790 to 1940 µm). Accordingly, the initial sizes of the cross-linking fronts are comparable as well and require accordingly almost the same time to reach the core of the spheres. The cross-linking fronts in replicates with 40 mM BaCl₂ need 120 to 150 seconds to reach the core and are also strongly dependent on the sphere’s diameter.

Overall, the tendency of cross-linking velocities of discs and spheres is comparable (Figure 4b and Figure 5d). The velocity of gelation in this work is defined as reduction of ungelled core and negative. The gelation results for spheres show the fastest cross-linking front for 40 mM BaCl₂ with a subsequent decrease with decreasing BaCl₂ concentration (40 mM > 20 mM > 10 mM, Figure 5d).

3.2. Numerical Solution

A time lapse sequence of the obtained simulation results for each cross-linker concentration (10 mM, 20 mM and 40 mM) is depicted in Figure 6, analogue to Figure 3. As previously mentioned, there is currently no shrinking behavior implemented in the model, leaving the alginate discs at the same diameter throughout the gelation process. These simulations are completed in similar amounts of time as the experimental tests, with 10 mM being slightly slower (3 min difference) and 20 mM and 40 mM being slightly faster (3 min and 4 min difference respectively). This difference, however, is mostly caused at the end of the simulation, when the gelation front shows an increase in speed, which does not happen in the experiments (see Figure 7). Otherwise, the two phases of fast front propagation followed by a slower linear phase are also present here.

The simulation results of gelation in Figure 7 are also available in the Supplemental Information as video sequences (Alginate-Disc_gelation_sim_10mM.avi, Alginate-Disc_gelation_sim_20mM.avi and Alginate-Disc_gelation_sim_40mM.avi).

4. Discussion

4.1. Alginate Gelation

In this study, the gelation kinetics were investigated using an image-based approach and a numerical model created based on the chemical reaction and the collected data. The experimental data reveal a correlation between the cross-linker concentration, the velocity of the traveling gelation front, and the overall gelation time. For the first time, the gelation kinetics of UHV-alginates were analyzed in detail using the described optical approach. In the literature, the general phenomenon is described for low-viscosity alginates [22], but no systematic study for UHV-alginates with a barium concentration lower than 50 mM had been conducted so far. In this study, the shape of the alginate was disc-like, as this allowed for higher throughput and better analysis. Although our preliminary data suggest that this kind of analysis is also feasible for sphere-like, free-floating alginates, it requires higher effort for analysis. The gelation kinetics reveal a strong dependency on the cross-linker concentration used. The main driver of this observation may be the diffusion of the barium ions through the alginate network. The greater the initial cross-linking concentration, the higher the gradient from the outside to the interior of the hydrogel. Due to the high concentration, barium is available in excess for the gelation reaction, and in consequence, no slowdown of the gelation front occurred in the linear phase. The approach described in this study is furthermore capable of analyzing the influence of small cargo (cells or organoids) or debris (protein agglomerates) on gelation behavior.

The data also suggest shrinkage of the alginate discs, probably caused by syneresis (the release of water during gelation). This shrinkage is not yet accounted for in the current numerical model. Since shrinkage (or syneresis, respectively) provides additional insights into the gelation kinetics, especially when two different alginate species are used. Alginates from LN contain more flexible heterogeneous MG/GM blocks, which enable coiling of polymers during gelation leading to shrinkage of the hydrogel. Syneresis in alginates is complex and the influence of the different components of gelation must be taken into account [55]. In further studies, the alginate composition, concentration of the cross-linking agent, and shrinkage will be studied in detail to further enhance the current numerical model.

The course of alginate gelation in spheres is strongly dependent of barium concentration. This observation was expected since a in a higher barium-concentrated solution, more ions are available to cross-link the alginate and the time, in which an ion finds a spot to cross-link a G-monomer is shorter. Therefore, it diffuses the fastest to the core of the UHV-alginate sphere. At the lowest concentrated BaCl2 solution, the availability of barium ions is lower, and it takes more time to find a cross-linking spot; therefore, it takes more time for the cross-linker to reach the core of the UHV-alginate sphere.

Overall, it seems that the optical observation of gelation kinetics is independent of the shape and volume of the alginates. Therefore, it will serve as a powerful and user-friendly experimental tool to further enhance the numerical model.

4.2. Layer-Wise Gelation

The microscopic images reveal a layer-like structure of the alginate discs. Over time, the number of layers at the outer surface increases. Similar structures have been described by Ehrhart et al. [56] and Jeong et al. [57]. In both studies, the depot of multivalent cations was the gelled alginate bead itself, which was suspended in liquid alginate. Thin alginate layers formed during gelation due to the diffusion of multivalent cations from the bead to the bead-alginate interface. A shell-wise gelation of alginate beads was also described in a review by Leong et al. [58] and Voo et al. [59]. In their study, Voo et al. found that highly concentrated (10% w/v) low-viscosity (LV) alginates form a multi-layer internal structure. Interestingly, this structure differs from low-concentrated (2% w/v) LV-alginates. Low-concentrated LV alginates resemble a more sponge-like structure with small cavities (pores). The different structures and arying degrees of homogeneity can be explained by the relative diffusion rate of alginate and the multivalent cation [59].

The experiments in this study were not designed to reveal the mechanism of layer formation during the gelation of UHV-alginates. However, initial interesting findings were made, indicating that such layer formation can be observed non-invasively during gelation with phase contrast microscopy (Figure 8a). In the given example, the thickness of the adjacent layers was approximately 15 to 30 µm for the 10 mM condition. These initial findings were summarized in Figure 8b. It can be hypothesized that the layer formation of UHV-alginates differs from the findings of LV-alginates, as polymers with higher molecular weight lead to differences in layer thickness.

For a detailed analysis and correlation with crosslinker concentration, alginate composition, and alginate concentration, a microscopic analysis with higher resolution and magnification must be conducted in further studies. By using alginates with different interactions with multivalent cations, these layers could help to understand the diffusion rates of the cations and polymers in UHV-alginates. In this context, future models will also consider the step-wise interaction of barium ions with alginate chains as presented for calcium by Zhou et al in a population balance model [60]. This model will be integrated and combined with the findings of layer-wise gelation to more precisely project the properties of the alginate gels at the hydrogel’s surface.

4.3. Numerical Simulation

Overall, the presented model achieves a satisfactory level of agreement with the experimental data and offers a strong starting point for subsequent improvement and expansion.

The sudden acceleration towards the end of Figure 7 could be due to variations in mesh refinement between the center and the outer rim of the circular domain. These adjustments may also help resolve the formation of steps observed in simulation results at smaller cross-linker concentrations. However, both issues are still being actively investigated. The spikes or abrupt increased at the very end of the graphs result from the automated front calculations detecting minor movements and simulation artifacts, as the primary reaction front has nearly dissipated.

Computation time can vary significantly depending on model parameters and the size, complexity and refinement of the used geometry mesh.

5. Conclusion and Outlook

In addition to the gelation of perfect hydrogels without any cells, cargo, or air bubbles, the experiments conducted also revealed the behavior of the cross-linking front after passing through some inclusions (debris). The acquired images revealed some swirling structures of the cross-linked hydrogels directly after the gelation front hits the object (data not shown). These structures, especially in high cell density immobilization applications such as bioprinting, must be further studied to understand the mechanical stability of the entire construct. Furthermore, the numerical model in this study, based on optical observations of the gelation and stoichiometry, serves as a solid foundation for a framework to optimize gelation properties. In the next step, the model will be refined by integrating nanoindentation data to project and predict the Young’s modulus of the hydrogels to be produced. By integrating this data in the numerical model, a prediction tool (e.g. as a mobile app) for the generation of tissue of interest can be implemented that will supervise cell biologists as well as material scientist to produce e.g. artificial cell niches for the applications in tissue engineering.