Bioink: Calcium Diffusion Modelling

This page details the computational approach for simulating diffusion and crosslinking processes in alginate-based ink formulations.

Introduction to the Hydrogel Crosslinking Mechanism

We were interested in computationally modeling the crosslinking of sodium alginate when exposed to calcium chloride within our bioink samples. This approach aims to quantify calcium ion diffusion based on a defined set of parameters, which will guide optimization of our bioink composition protocols. This will help ensure structural integrity of our printed structures.

Hydrogel crosslinking within our bio-ink is visualized in Figure 1 below:

Calcium ions (Ca2+) from the calcium chloride solution diffuse into the gel, where some are captured by sodium alginate through crosslinking, forming a solidified gel zone. As more ions continue to diffuse, the already gelled region creates a physical barrier, slowing diffusion. This interaction between Ca2+ ion movement and gelation creates a crosslinking zone. Additionally, sand present in the bioink blocks calcium ion diffusion and crosslinking networks.

ASSUMPTION 1: Diffusion Modelling with Alginate-Regolith Bio-Ink Only

We focused on diffusion modeling within alginate-regolith ink, as these are the primary components of our bioink. To simplify the problem, we excluded the effects of CMC, UTEX2973, and its growth media on diffusion and crosslinking.

DBTL 1 - Microscope-based Diffusion + Crosslinking Modelling

Our Approach

Our initial modeling approach focused on observing in vivo cross-linking in ink samples to develop intuition about the underlying processes. We conducted microscopy-based experiments 2025.06.11 Calcium Diffusion Under Microscope (1) and 2025.06.12 Calcium Diffusion Under Microscope (2) to track calcium ion diffusion, drawing from [1]. This methodology was selected for its accessible computational strategy and use of microscopy to visualize regolith particle size and distribution.

Computational Strategy

The model shows how the crosslinking agent diffuses in three dimensions by calculating small changes over time and distance. They specifically do this by using Equation 1, which demonstrates the dependence of crosslinking zone displacement $x_c(t)$ on the square root of the time elapsed.

$$x_c(t) = k\sqrt{t}$$

Equation 1 is obtained from the analytical solution to the diffusion equation in Cartesian coordinates. The authors also demonstrate the linear dependance of $x_c$ on $\sqrt{t}$ for solutions in spherical and cylindrical coordinates approximated with a second-order Taylor series expansion. The consistency across coordinates shows that the square root relationship between crosslinking distance and time holds true in different shapes, making the model more reliable.

$$k = \sqrt{\frac{4C_0 D}{C_g}}$$

$C_o$: initial concentration of Ca2+

$D$: diffusion coefficient of the free Ca2+ ions in the cross-linked scaffold

$C_g$: concentration of binding sites of alginates that can react with Ca2+ ions

The constant $k$, whose analytical basis is established from the diffusion equation, can be obtained experimentally from the gradient of the line produced by plotting $x_c$ against $\sqrt{t}$ . This constant can then be used to compute the instantaneous diffusion time at a given penetration depth.

Methodology

The methodology aims to experimentally determine the value of $k$ for our bio-ink samples. This procedure primarily relies on the fact that cross-linked alginate turns opaque compared to its free form. The cross-linking front displacement increments $x_c(t)$ can then be measured by pixel analysis. A microscope is used to sample the cross-linking displacement increments, since the relationship in Equation (1) only holds for small penetrations and small times.

Results

ASSUMPTION 2: “Simple Diffusion with a Barrier” Modelling Approach

Based on Figure 3, we assumed calcium ion diffusion in the bio-ink follows simple diffusion with sand particles acting as barriers, since their dispersion was insufficient to classify the ink as porous media.

We tested two bioink samples: one with pure alginate and one with 30% earth sand. While the alginate-only sample showed clear crosslinking behavior, the sand-containing sample lacked observable opacity changes, preventing us from tracking the crosslinking front. Additionally, the crosslinker failed to diffuse beneath the microscope slide in the sand-alginate sample, leaving it soft and uncrosslinked.

Next Steps

After encountering failure during Testing, we experimented with increasing the crosslinking time and reducing the bio-ink volume applied to the slide from 200μL to 50μL. Despite these adjustments, persistent issues with ink viscosity and visual interference from sand particles led us to adopt an alternative modeling approach.

DBTL 2 - Volumetric-based Diffusion + Crosslinking Modelling

Our Approach

We modelled sodium alginate crosslinking using the approach from [2], which predicts calcium chloride absorption over time. This was ideal for our bioink containing sand, which blocks visual tracking. The model combines partial and ordinary differential equations to describe ion diffusion and gelation in one dimension. After validating with a control experiment, we calibrated it to our sand alginate bioink using feedback from Prof. Michele Marino, enabling diffusion estimates across different extrusion geometries.

Computational Model

The computational model by Hajikhani et al. models the geometry of the problem as one-dimensional ([2]). This means that the calcium ions diffuse along a line of length $L$, seen in Figure 1. A variable $x$ is introduced, which represents the position along this line, i.e. $x \in [0,L]$. The variable $t$ is also introduced, which is the time elapsed, $t \in [0,t_{end}]$.

The authors then define two functions of $x$ and $t$, which are the calcium ion concentration and gelation degree.

Calcium ion concentration: $$c = c(x,t)$$ Gelation degree: $$\alpha = \alpha(x,t)$$

The gelation degree is mathematically represented by the ratio of the concentration of solidified gel ($c_g$) to the initial concentration of free alginate ($c_a$). This ratio serves as a quantitative indicator of the extent of alginate solidification.

$$\alpha = \frac{c_g}{c_a} \in [0,1]$$

When $\alpha = 0$, the concentration of solidified gel is zero ($c_g = 0$), indicating no cross-linking has occurred and the gel remains liquid. When $\alpha = 1$, the concentration of solidified gel equals the initial concentration of free alginate ($c_g = c_a$), meaning all alginate chains are fully crosslinked and the gel is completely solidified.

Diffusion Coefficient as a function of Gelation Degree ($\alpha$):

As shown earlier in Figure 1, the progression of the crosslinking zone inhibits Ca2+ diffusion. One way to model this is to parametrize the diffusion coefficient as a function of the gelation degree ($\alpha$). In other words, the physical state of the gel affects the diffusivity of calcium ions by the relationship below.

$$D(\alpha) = D_0 + (D_1 - D_0)\frac{e^\frac{-n\alpha} {\alpha_{gel}} - 1{}}{e^\frac{-n}{\alpha_{gel}} - 1}$$

Where:

• $D_0$: Diffusion coefficient of Ca2+ ions across the liquid gel ($\alpha = 0$).
• $D_1$: Diffusion coefficient of Ca2+ ion across the fully solid gel ($\alpha = 1$).
• $\alpha_{gel}$: The characteristic gel point defined by $0 \leq \alpha_g \leq 1$. It is a value of the gelation degree for the alginate to just start behaving as a solid, without being fully crosslinked.
• $n$ is a parameter which sets the rate of diffusivity change.

In our iHP meeting with Prof. Marino, we learned that the diffusion coefficient function $D(\alpha)$ is a postulation by the authors, which they found led to numerical results aligned closely with experiments.

Reaction-Diffusion Equation:

The reaction-diffusion partial differential equation in 1D is derived from a simplified form of the Mikkelsen-Elgsaeter model ([3]). An assumption made here is the diffusion of free alginate molecules out of the gel is negligible.

$$\frac{\partial c}{\partial t} - \frac{\partial }{\partial x}\left(D(\alpha)\frac{\partial c}{\partial x}\right) = -r =\underbrace{-N_c c_a\frac{\partial \alpha}{\partial t}}_\textrm{sink term}$$

Where:

• $N_c$: The stoichiometric coefficient for the number of Ca2+ ions used up per alginate-alginate crosslink The effects of crosslinking are accounted for with a “sink” term. Alginate cross-linking acts as a sink to free diffusing Ca2+ ions by depleting their concentration when within its lattice.

1. Gelation Reaction Equation: $$\frac{\partial \alpha}{\partial t} = K \frac{c}{c_a}(1 - \alpha)$$

The partial derivative on the left hand side of the equation is the evolution of the gelation degree which depends on:

• $c(x,t)$: the concentration of Ca2+ ions
• $c_a$: the initial concentration of free alginate
• $K$: an experimentally fitted model parameter associated with the reaction rate.

Computational Strategy

To solve the system of equations for our experimental conditions we used FEniCS, an open source finite element library PDE solver. Our full implementation is available on GitLab, but we have gone through the steps of the program.

Boundary and Initial Conditions:

The first step to our computational strategy was implementing the provided boundary and initial conditions. This was done using the FEniCS boundary and initial condition set-up for finite element analysis.

$$B.C. \begin{cases} c=\^{c} &\text{for } x = 0\\ \frac{\partial c}{\partial x} = 0 &\text{for } x = L \end{cases}$$ $$I.C. \begin{cases} c(x,0) = 0 &\text{for } x \neq 0 \end{cases}$$

Where:

• $\^{c}$ : The concentration of Ca2+, where $\^{c} = w_{ca}\^c_{ca}$
• $w_{ca}$: The weight fraction of Ca2+ in CaCl2 solution.
• $\^c_{ca}$: The initial concentration of Ca2+ at the boundary.

Discretization of the Temporal and Spatial Domain:

The time interval $t \in [0,t_{end}]$ is discretized into $n$ time intervals such that $\Delta t = t_{n} - t_{n -1}$. We iterate over $n$ to cover the total time. The spatial domain is $x \in [0,L]$, represented as a mesh, is split into $N$ nodes to create the finite element grid space.

1. Euler’s Method to solve Gelation Reaction Equation: The third step was to use Euler’s method to compute the next iteration $\alpha_n$ recursively using the previous value, $\alpha_{n-1}$. We use the general form of Equation 11 in our program, however our indexing starts from $n = 0$, with an assigned initial condition of $\alpha_0 = 0$. These gelation degree values are then used to compute the diffusion coefficient $D (\alpha)$

$$\alpha_n = \alpha_{n-1} + K\frac{c_{n-1}}{c_a} \left(1-\alpha_{n-1}\right)\Delta t$$

The values $\alpha_n$ and $\alpha_{n-1}$ are used for a finite difference approximation of the time derivative of the gelation degree, shown in Equation 12. This is used to assign the RHS sink term of Equation 7, when solving for the variational problem in the next step.

$$\frac{\partial \alpha}{\partial t} \approx \frac{\alpha_n - \alpha_{n-1}}{\Delta t}$$ FEniCS Finite Element Solver:

FEniCS uses the finite element method, which works with PDEs in their weak form. We need a trial function $u$, which is a symbolic placeholder for the unknown calcium concentration $c_n$ , and a test function $v$ to convert our PDE to this form.

We approximate the time derivative of the calcium ion concentration as:

$$\frac{\partial c}{\partial t} \approx \frac{c_{n} - c_{n-1}}{\Delta t}$$

Then substitute it into Equation 7, multiply by our test function $v$ and integrate over the domain to convert the PDE to its weak form. The vector notation is used for simplicity in FEniCS programming, even though our problem is in 1D and does not require it.

$$\int_{\Omega} \frac{c_n - c_{n-1}}{\Delta t}\, v \, dx-\int_{\Omega} \nabla\!\cdot\!\bigl(D_{\alpha}\,\nabla c_n\bigr)\, v \, dx= - \int_{\Omega} r_n\, v \, dx$$

The input to the solver function in FEniCS requires the weak form in Equation 14 to be expressed in it’s discrete variational form, shown in Equation 15.

$$a(u,v) == l(v)$$

After integrating by parts, rearranging, and replacing $c_n $ with our trial function $u$, Equation 14 can be written in the form of Equation 15 by:

$$\underbrace{\int_{\Omega} \frac{1}{\Delta t}\, u\, v \, dx\;+\; \int_{\Omega} D_{\alpha}\, \nabla u \cdot \nabla v \, dx}_{a(u,v)}\;=\;\underbrace{\int_{\Omega} \frac{1}{\Delta t}\, c_{n-1}\, v \, dx\;-\; \int_{\Omega} r_n\, v \, dx}_{l(v)}$$ Computing volume at the boundary:

The last step of the program is to compute the volume of calcium chloride absorbed through the top surface of the gel over the time elapsed.

$$V_c(t)=\frac{A}{w_{ca} \rho_{\mathrm{CaCl}_2}}\int_{0}^{t} D\,\left.\frac{\partial c_c}{\partial x}\right|_{x=0}\, dt$$

Where:

• $A$: The cross-sectional area of the vial. Here it is a circle, so $A = \pi r^2$, where r is the radius of the vial.
• $\rho_{\mathrm{ca}}$: The solution density of calcium chloride.

Methodology

The experimental protocol can be found in 2025.06.21 Calcium Diffusion in Vials (Earth Sand) and 2025.09.21 Calcium Diffusion in Vials (MGS-1). We simply crosslinked 1 mL of Earth Sand or MGS-1 alginate gel with 1 mL of 100 mM calcium chloride over different times as shown below:

The difference in weight before and after crosslinking was calculated and converted to volume of calcium chloride absorbed through the following equation:

$$\text{v}_{CaCl2} =\frac{w_{f} - w_{i}}{\rho_{CaCl_2}}$$

Where the difference in weight is divided by density of calcium chloride. We approximated the density as 1.007 g/mL. The assumption made is that the aqueous calcium chloride solution is uniformly mixed, so the fraction of calcium chloride molecules per volume of crosslinker absorbed by the gel can be considered constant.

Model Validation

After executing our computational strategy, we wanted to validate our results to ensure the model was working properly. To do this, we fitted model parameters from [2], in order to recreate Figure 3(a) from their paper. The computational model is run for two cases of $D(\alpha)$.

Case where $D(\alpha)= D_0$:

Assumes constant diffusion coefficient $D_0$. Keeping the diffusion coefficient constant is a much easier starting point for validating results than the original expression in Equation 6. The downside is that it does not fully capture the effects of cross-linking by neglecting the changes in diffusion coefficient from gelation. This results in a curve that deviates away shown in Figure 6.

Case where $D(\alpha)$ takes the form in Equation 6:

The use of this expression in our computation results in a volume output that is much closer to experimental values.

The full list of parameters used for the above can be found in our github repository. Comparing our results to [2] Figure 3(a) reveal close alignment in the absorbed volume change profiles for both $D_0$ and $D_{\alpha}$. Although some differences in values are observed, such as near the boundaries, these may be due to how the solvers handle discretization or error propagation. They may also be attributed to mesh or temporal resolution variations. Still, the strong agreement between our results and those referenced indicate our approach was successful at implementing their model.

Next, we wanted to validate the model with our own results. To do this, we ran our methodology for our alginate-only (0% sand) bioink samples which was done on 2025.08.04 Calcium Diffusion in Vials (Alginate Control Repeat).

The full list of parameters used for the above can be found in our github repository. After determining the absorbed volume of calcium chloride from our measurements, we then overlayed the data points onto our numerical results shown in Figure 8. The model parameters are adjusted for the geometry of our vials, as well as differences in reagent concentration. The fitted $D(\alpha)$ model had $R^2 = 0.84$, demonstrating good agreement with the experimental data. However, we found our $RMSE = 29.2 \%$, which most likely is attributed to random error from our experiments, or model deviation.

Results

We repeated our experimental methodology for the 30% sand and 20% MGS + 3.5% CMC bio-ink samples. This decision was informed by our iHP meeting with Prof. Marino, who recommended adjusting parameters in Equation 6 to better align our model with the experimental data. The impact of regolith incorporation on the crosslinking behavior of alginate-based bioinks remains largely unexplored due to the novelty of the approach. However, we noted that sand particles act as barriers to Ca2+ ion diffusion, likely influencing the parameters $D_0$, $D_1$ and $\alpha_{gel}$, which characterize the bio-ink’s solidification dynamics.

The full list of parameters used for the above can be found in our github repository. For the 30 wt% earth sand and 9 wt% CMC alginate gel measurement, we attained good results by fixing $D_0$ and $\alpha_{gel}$, while setting $D_1 = 0.25D_0$. Our calibrated model $D(\alpha)$ curve had an $R^2 = 0.97$, indicating a strong trend fit between the model and experimental data. Additionally, the $RMSE = 7.8$%, suggesting that the model’s predictions closely matched observed values with relatively low error.

Parameter fitting for our 20 wt% MGS-1 and 3.5% wt% CMC alginate gel was more challenging. We found that fixing $D_0$, while setting $D_1 = 0.05D_0$ and $\alpha_{gel} = 0.8$, yielded the best calibrated model fit, with an $R^2 = 0.47$ and $RMSE = 16.49\%$. Unfortunately, parameter fitting was done heuristically rather than through formal optimization. However, our choices were guided by physical reasoning: the diffusion coefficients were lowered to reflect reduced Ca2+ diffusion from particle obstruction, and $\alpha_{gel}$ was increased to account for the higher ion concentration required to initiate gelation in the presence of these barriers.

Future Directions

For future modelling, we would like to study how embedded particles quantitatively affect Ca2+ diffusion, or cause lattice destabilization. This would inform our calibration approach to $D_0$, $D_1$, and $\alpha_{gel}$ A potential future step might be to train a computer model to estimate optimal parameter values. The calibrated model will then be tested on extruded ink samples and hopefully extended to two and three dimensions to explore living building material geometries that enhance crosslinking and MICP, ensuring structural integrity.