(C) PLOS One
This story was originally published by PLOS One and is unaltered.
. . . . . . . . . .



A computational algorithm for optimal design of a bioartificial organ scaffold architecture [1]

['Martina Bukač', 'Department Of Applied', 'Computational Mathematics', 'Statistics', 'University Of Notre Dame', 'South Bend', 'Indiana', 'United States Of America', 'Sunčica Čanić', 'Department Of Mathematics']

Date: 2024-12

Abstract We develop a computational algorithm based on a diffuse interface approach to study the design of bioartificial organ scaffold architectures. These scaffolds, composed of poroelastic hydrogels housing transplanted cells, are linked to the patient’s blood circulation via an anastomosis graft. Before entering the scaffold, the blood flow passes through a filter, and the resulting filtered blood plasma transports oxygen and nutrients to sustain the viability of transplanted cells over the long term. A key issue in maintaining cell viability is the design of ultrafiltrate channels within the hydrogel scaffold to facilitate advection-enhanced oxygen supply ensuring oxygen levels remain above a critical threshold to prevent hypoxia. In this manuscript, we develop a computational algorithm to analyze the plasma flow and oxygen concentration within hydrogels featuring various channel geometries. Our objective is to identify the optimal hydrogel channel architecture that sustains oxygen concentration throughout the scaffold above the critical hypoxic threshold. The computational algorithm we introduce here employs a diffuse interface approach to solve a multi-physics problem. The corresponding model couples the time-dependent Stokes equations, governing blood plasma flow through the channel network, with the time-dependent Biot equations, characterizing Darcy velocity, pressure, and displacement within the poroelastic hydrogel containing the transplanted cells. Subsequently, the calculated plasma velocity is utilized to determine oxygen concentration within the scaffold using a diffuse interface advection-reaction-diffusion model. Our investigation yields a scaffold architecture featuring a hexagonal network geometry that meets the desired oxygen concentration criteria. Unlike classical sharp interface approaches, the diffuse interface approach we employ is particularly adept at addressing problems with intricate interface geometries, such as those encountered in bioartificial organ scaffold design. This study is significant because recent developments in hydrogel fabrication make it now possible to control hydrogel rheology and utilize computational results to generate optimized scaffold architectures.

Author summary Bioartificial organ is an engineered device made of living cells and a biocompatible scaffold that can be implanted or integrated into a human body to replace a natural organ, or to augment a specific organ function. The biocompatible scaffold serves as the structural foundation, embedding cells and growth factors to form substitute tissue. A key challenge in bioartificial organ design is ensuring the long-term survival of transplanted cells. Specifically, a critical challenge is creating a hydrogel architecture with ultrafiltrate channels that act as pathways to deliver oxygen and nutrients to the cells, supporting their sustained viability. In this work, we developed a mathematical and computational framework to investigate optimal hydrogel architecture and oxygen supply to transplanted cells. Our findings show that a hexagonal channel architecture significantly outperforms both the commonly used vertical channels and branching channel designs, providing a substantial increase in oxygen delivery to the cells. In the hexagonal architecture, oxygen concentration remains above the critical threshold for hypoxia throughout the scaffold, which is not achieved with other designs. This is particularly important given recent advances in hydrogel fabrication, allowing for precise control of hydrogel elasticity and rheology, making these results valuable for developing improved scaffold designs.

Citation: Bukač M, Čanić S, Muha B, Wang Y (2024) A computational algorithm for optimal design of a bioartificial organ scaffold architecture. PLoS Comput Biol 20(11): e1012079. https://doi.org/10.1371/journal.pcbi.1012079 Editor: Jason M. Haugh, North Carolina State University, UNITED STATES OF AMERICA Received: April 14, 2024; Accepted: October 16, 2024; Published: November 11, 2024 Copyright: © 2024 Bukač et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: The computational codes are publicly available on GitHub at the following link: https://github.com/qcutexu/StokesBiotSolver. Funding: This work was supported by the National Science Foundation (DMS-2208219 to MB; DMS-2205695 to MB; DMS-224700 to SC; DMS-2011319 to SC; DMS-2247001 to YW). This work was also supported by Hrvatska Zaklada za Znanost (IP-2022-10-2962 to BM) and Ministry of Science Education and Youth of Croatia (Croatia-US Bilateral Grant “The mathematical framework for the diffuse interface method applied to coupled problems in fluid dynamics” to BM). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

1 Introduction Bioartificial organ is an engineered device made of living cells and a biocompatible scaffold that can be implanted or integrated into a human body–interfacing with living tissue–to replace a natural organ, or to duplicate or augment a specific organ function [1]. Biocompatible scaffold is a biocompatible base material in which cells and growth factors are embedded to construct a substitute tissue, which can be used in e.g. bioartificial organ design [2]. For example, a biocompatible scaffold used in the bioartificial pancreas design is a biocompatible hydrogel (e.g. agarose gel) which houses transplanted pancreatic islets (conglomerations of pancreatic cells) that contain the insulin-producing β-cells. Hydrogels are hydrophilic polymers that are poroelastic, and that have the ability to absorb a large volume of fluid, which makes them particularly suitable materials for biomedical applications. In bioartificial pancreas design presented in [3], the biocompatible hydrogel containing the insulin-producing cells is encapsulated in a nano-pore, semi-permeable membrane (filter), see Fig 1, and the resulting bioartificial organ is connected to the host’s cardiovascular system via an anastomosis graft (a tube) for the advection-enhanced nutrients delivery to the cells, and insulin distribution away from the cells. PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 1. Schematic illustration of immunoprotective encapsulation device containing human embryonic stem cell-derived β clusters (hES-βC). https://doi.org/10.1371/journal.pcbi.1012079.g001 A major challenge in the encapsulation-based bioartificial organ design is maintaining the survival of transplanted cells within the bioartifical organ scaffold for an extended period of time. In particular, long-term viability of pancreatic cells is directly affected by the sufficient access to nutrients for survival, of which oxygen is the limiting factor [4]. One of the goals of this manuscript is to develop the mathematical and computational approaches to investigate oxygen supply to the transplanted cells by studying hydrogel architecture and optimal design of ultrafiltrate channels within the encapsulated poroelastic hydrogels for advection-enhanced reaction-diffusion processes of oxygen transport to the cells. This is significant because recent developments in hydrogel fabrication make it now possible to control hydrogel elasticity and hydrogel rheology using, e.g., the approaches employed and reviewed in [2, 5]. In this manuscript we study oxygen-carrying blood plasma flow and oxygen concentration in three different scaffold architectures, whose geometries are inspired by different biological flow-nutrients scenarios. The first geometry consists of vertical ultrafiltrate channels drilled through a hydrogel, which is the simplest, and a standard procedure used in the design of bioartificial pancreas [5–7]. The second geometry consists of the branching channels, inspired by the architecture of the branching vessels in the human body. The third geometry consists of the channels surrounding hexagonal pockets, which we refer to as the hexagonal geometry. This was inspired by the biological (epithelial) tissues in which interstitial fluid flows through a network of irregularly arranged interstices between hexagonally shaped cells, which supports their structural and functional integrity [8]. To work with comparable fluid flow scenarios, the three geometries are generated so that the total fluid channels’ volume (i.e., the area of the 2D channels) is the same in all three geometries. We are interested in a scaffold architecture with the geometry that provides concentration of oxygen that is as uniform throughout the scaffold as possible and above a minimum concentration below which hypoxia occurs. For this purpose we introduce two sets of models, one for the fluid flow and one for oxygen concentration in the scaffold. The scaffold is modeled as a 2D domain containing a network of channels separating the regions of poroelastic medium describing a hydrogel. See Fig 2 below. In the last part of the manuscript we introduce a 3D computational model to show that for the “optimal” geometry, the 2D simulations capture the main features of flow and oxygen concetration. The flow of blood plasma through the network of channels is modeled by the time dependent Stokes equations for an incompressible, viscous fluid. The channel flow is coupled to the filtrating flow in the poroelastic hydrogel, which is modeled by the Biot equations [9]. Two sets of coupling conditions at the interface between the channel flow and hydrogel poroelastic medium are employed: the kinematic and dynamic coupling conditions. The kinematic coupling conditions describe continuity of normal (vertical to the interface) velocity and a slip in the tangential component of the velocity, known as the Beaver-Joseph-Saffman condition. The dynamic coupling conditions state the balance of forces at the interface. Linear coupling is considered, i.e., the coupling conditions are evaluated at a reference interface. Subsequently, the resulting fluid velocity is used as an input data in an advection-reaction-diffusion model for oxygen concentration within the scaffold. PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 2. An example of a domain Ω = Ω F ∪ Ω P demonstrating the two inlets An example of a domain Ω = Ω∪ Ωdemonstrating the two inletsat the top and the two outletsat the bottom in 2D (left) and 3D (right). https://doi.org/10.1371/journal.pcbi.1012079.g002 To solve the coupled problems computationally, we implemented a diffuse interface approach [10–12]. The diffuse interface approach has several advantages over the classical sharp interface approach. In the diffuse interface approach, the mesh nodes do not have to be aligned with the interface, whose location is now captured using a phase-field function. Instead, the problem unknowns are defined on the entire domain, and the mesh remains fixed. This is particularly useful in cases where the interface between the two regions is difficult to determine exactly (e.g. “fuzzy” boundaries between regions obtained from medical images), or when the geometry of the interface is complex, leading to inaccurate approximations of the boundary integrals. Furthermore, advantages of this approach are apparent in problems with moving geometries and in the development of geometric optimization algorithms for optimal channel architecture design, which we plan to do next. In these cases a diffuse interface approach avoids the creation of complicated meshes following optimizing scaffold geometries, which would be the case in the classical sharp interface approaches. In this work we develop a diffuse interface approach to solve both the fluid and oxygen problems. In both cases, a finite element method is used for spatial discretization, and a backward Euler method for time stepping. We perform a verification of our numerical diffuse interface approach by comparing the solutions of the diffuse interface problem with the already validated solutions of a sharp interface approach, showing excellent agreement. Additionally, to justify the computationally less expensive 2D simulations, for the optimal geometry we perform the full 3D simulations and show that the 2D simulations we used to find the optimal geometry, capture all the important features of plasma flow and oxygen concentration in the corresponding 3D scaffold. The diffuse interface approach is then used to investigate flow, pressure, and oxygen concentration in the three main geometries, discussed above. As mentioned earlier, the diffuse interface approach developed here is particularly suitable from the computational standpoint to study a number of different interface geometries. We find that the ultrafiltrate channel network with the hexagonal geometry provides the best solution in terms of oxygen concentration distribution and magnitude. Oxygen concentration in the hexagonal case is most uniformly distributed throughout the entire scaffold, and the values of concentration are all well above the critical value of oxygen below which hypoxia occurs. This is not the case for the other two geometries considered here (vertical channels and branching channels). We then further investigate the reasons why the hexagonal geometry performs best, and conclude that the relative angle of the channels with respect to the flow direction is one of the main contributing factors to higher Darcy velocity within the poroelastic hydrogel, and consequently the higher concentration values in between the ultrafiltrate channels. The results obtained in this manuscript can be used in optimal scaffold design by implementing hydrogel manufacturing techniques recently developed in [2, 5]. This manuscript is organized as follows. In Section 2 we introduce the model equations in strong/differential and weak/integral formulations. In Section 3 the diffuse interface method is discussed, and the fluid and oxygen problems are presented in the weak (integral) diffuse interface formulations. Numerical discretization based on a finite element approach is also discussed in this section. The computational setting, which includes the computational geometries described by the phase-field function, and parameter identification, is presented in Section 4. In this section we also show the verification of the numerical solver by comparing the diffuse interface solution with a monolithic solver utilizing a sharp-interface approach. Numerical results are presented in Section 5. In Section 6, we present 3D simulations of the optimal geometry demonstrating strong agreement with the 2D simulations. Conclusions are presented in Section 7.

2 Methods I: The mathematical model To model the entire scaffold, we introduce a bounded domain , with d = 2, 3, consisting of a fluid region, Ω F , corresponding to the blood or blood plasma flow region bringing oxygen and nutrients to the transplanted cells, and a poroelastic solid region, Ω P , corresponding to the poroelastic medium such as a biocompatible hydrogel, which hosts the transplanted cells. The entire bioartificial organ scaffold is then comprised of the two regions so that and Ω F ∩ Ω P = ∅ (see Fig 2). The fluid flowing through the channels comprising the region Ω F enters the poroelastic medium Ω P across the interface Γ, which is the interface between Ω F and Ω P , defined by . Fig 2 shows a sketch of a 2D and a 3D domain in which the fluid domain Ω F consists of two horizontal channels corresponding to the “inlet” and “outlet” channels, namely, they contain the inlet boundary and the outlet , respectively. The two horizontal channels are connected via three vertical channels describing the simplest scaffold architecture which has been used in several research studies on bioartificial organ design [5–7]. The 3D geometry is obtained by extruding the 2D geometry. The vertical ultrafiltrate channels are drilled to enhance advection-dominated oxygen and nutrients supply to the cells residing in Ω P . This is just one, and the simplest example of the scaffold architecture that we will consider in this manuscript. To model the interaction between the flow of blood plasma in Ω F and the filtration of blood plasma through Ω P , we will use a fluid-poroelastic structure interaction model, which we present next. We assume that the fluid-poroelastic structure interaction problem is linear, and that the displacement is small enough so that the domain remains fixed. 2.1 Fluid-poroelastic structure interaction We assume that Ω F contains a viscous, incompressible, Newtonian fluid, such as blood plasma, whose dynamics can be described by the Stokes equations: (1a) (1b) where v stands for the fluid velocity, ρ F is density, and σ F is the Cauchy stress tensor defined by: (2) where p F represents the fluid pressure, μ F is the fluid viscosity, and D(v) = (∇v + ∇vτ)/2 is the strain rate tensor. Namely, Eq (2) is a constitutive law describing a Newtonian fluid such as blood plasma. Velocity v will be used in the advection-diffusion-reaction problem to convect the oxygen concentration in Ω F . The poroelastic medium, such as agarose gel used in the design of a bioartificial pancreas, see [4], can be described by the well-known Biot model, which has been used to describe hydrogels in [13]. The Biot model consists of the following equations: (3a) (3b) (3c) Eq (3a) describes the elastodynamics of the elastic skeleton, i.e., the solid phase in the Biot model, and is given in terms of the displacement of the elastic skeleton, denoted by η, from its reference configuration Ω P . Constant ρ S denotes density of the poroelastic matrix, and γ E is a spring coefficient. For the simulations performed in 2D, adding the term γ E η accounts for the three-dimensional elastic energy effects, and helps to keep the regions that are not connected from drifting away. More precisely, this term helps to prevent spurious numerical solutions where sections of Ω P might be excessively pulled in the direction of vertical fluid flow. Such unrealistic large vertical displacements do not occur in 3D models due to the presence of an elastic restoring force in the third spatial dimension. In 3D simulations this term is not needed, so in that case we set γ E = 0. The total Cauchy stress tensor σ P for the poroelastic region is defined by where α is the Biot-Willis parameter accounting for the coupling strength between the fluid and solid phases and σ E denotes the elastic stress tensor, described by the Saint Venant-Kirchhoff constitutive model as: where λ E and μ E are Lamé parameters. Eq (3b) is the Darcy law, where p P is the fluid pore pressure, q is the Darcy velocity, and κ is the permeability. Eq (3c) is the storage equation for the fluid mass conservation in the pores of the matrix. It describes the coupled behavior of fluid flow and solid deformation in a porous medium. In this context, the quantity c 0 p P + α∇ ⋅ η represents the total fluid content, where c 0 is a constant related to the fluid’s compressibility. The mass conservation equation incorporates the effects of both fluid and solid compressibility to ensure that the total mass is conserved, considering how the fluid pressure changes impact the solid matrix and vice versa. The Darcy velocity q will be used in the advection-diffusion-reaction problem to convect the oxygen concentration in Ω P . Coupling conditions. To couple the fluid flow model with the Biot poroelastic medium model we use a set of four coupling conditions: two kinematic coupling conditions and two dynamic coupling conditions, which are evaluated at the location Γ of the interface between the two models. To state the coupling conditions introduce n F to denote the outward unit normal vector to ∂Ω F . At the interface Γ the following coupling conditions hold: (4a) (4b) (4c) (4d) where for a vector v, we define a projection onto the local tangent plane on Γ as: The first two conditions provide information about kinematic quantities such as velocities (kinematic coupling conditions), and the second two conditions provide information about coupling of stresses/forces (dynamic coupling conditions). More precisely, Eq (4a) describes the fluid mass conservation across the interface, i.e., the normal components of the free fluid velocity v relative to the velocity of motion of the poroelastic matrix ∂ t η is equal to the Darcy velocity q across the interface. Eq (4b) describes the Beavers-Joseph-Saffman condition (4b) with slip rate γ > 0, namely, the tangential component of the free fluid velocity slips at the interface with the slip rate proportional to the fluid shear stress (σ F n F ) τ . Eqs (4c) and (4d) describe the continuity of pressures and total stresses at the interface. Boundary and initial conditions. Problem (1a)–(4d) is supplemented with boundary and initial conditions. In our notation, we represent the exterior boundary as Γext = ∂Ω. The exterior boundary is divided into three distinct parts: the fluid inflow boundary, the fluid outflow boundary, and the impenetrable part of the boundary, . For the fluid problem, we impose the following boundary conditions: where v in is a prescribed velocity, specified in Section 4. For the poroelastic structure, we impose no flow and zero displacement at the external boundaries: (5) (6) Initially, we assume that the fluid is at rest and that the deformable poroelastic structure is in its reference configuration. Thus, we have: (7) Weak formulation. To define the weak formulation of problem (1a)–(7) we introduce the following function spaces. Given an open set S, we consider the usual Sobolev spaces Hk(S), with k ≥ 0, and introduce the following function spaces: Here corresponds to the test space for the fluid velocity, is associated with the solution space for the fluid velocity, and is associated with the function space for the poroelastic matrix displacement. The spaces H1(Ω P ) and L2(Ω F ) are associated with the solutions spaces for the fluid pressures in Ω p and Ω F respectively. We say that is a weak solution if for every , the following equality is satisfied in : (8) Here, denotes the space of distributions on (0, T), which is the dual space of the space of test functions . We note that the weak formulation is obtained using the primal formulation for the Biot problem, i.e., Eqs (3b) and (3c) have been combined so that Eq (3c) is written only in terms of p P and η. The Darcy velocity can be computed by postprocessing using (3b). Once the fluid velocity v and Darcy velocity q are computed from the fluid-poroelastic structure interaction problem specified above, we can use this velocity information to formulate an advection-reaction-diffusion problem describing oxygen transport in the bioartificial organ consisting of ultrafiltrate channels Ω F and the poroelastic medium Ω P containing the cells. 2.2 Advection-reaction-diffusion To model the transport of oxygen, we use an advection-reaction-diffusion model describing oxygen transport in the fluid domain Ω F and in the poroelastic medium Ω P . Oxygen transport in the human vascular system and tissues has been studied by many authors [14–18], and we adopt the approach from [15] to study oxygen transport in the scaffold: (9) where c is concentration of oxygen, and D is a diffusion coefficient equal to D F in Ω F and to D P in Ω P . The velocity u is equal to v in Ω F , and to q in Ω P . The final time is T t , which we note is different from the final time T used in the fluid-poroelastic structure interaction problem. The reaction term on the right-hand side is active in Ω P only and it accounts for the consumption of oxygen in the scaffold. In particular, this reaction term depends on the maximum oxygen consumption rate, R max < 0, the Michaelis-Menten constant, c MM , corresponding to the oxygen concentration when the consumption rate drops to 50% of its maximum [15], and the function , defined by accounts for the regions of the tissue/poroelastic matrix where the oxygen concentration falls below a critical concentration c CR below which necrosis is assumed to occur [15, 19]. The values of the parameters are all specified in [15, 19]. Boundary and initial conditions. Eq (9) is supplemented with the following boundary conditions: (10a) (10b) (10c) where c in is a given quantity, specified in Section 4. Thus, these boundary conditions say that we have a prescribed oxygen concentration at the inlet, zero oxygen concentration at the top boundary of the scaffold between the inlet regions, and zero diffusive oxygen flux at the outlet. Initially, the concentration is set to zero: (11) Weak formulation. To write the weak formulation of problem (9)–(11) we introduce the following function spaces: Here is associated with the test space for oxygen concentration and is associated with the solution space for c. We say that is a weak solution if c ≥ 0 and if for every , the following equality is satisfied in : (12) We have now specified two problems: a fluid-poroelastic structure interaction problem and an advection-reaction-diffusion problem, that we would like to solve for the fluid velocity and poroelastic structure displacement, and for oxygen concentration. The plan for this manuscript is to investigate three different scaffold architectures, motivated by biological structures, and numerically test which one provides the scaffold architecture with oxygen concentration that is closest to the uniform distribution of oxygen and is above the known minimal value c opt for which uninhibited insulin production by the β-cells is guaranteed [20]. See Section 5. For this purpose we developed two numerical methods, one for the fluid-poroelastic structure interaction problem, and one for the advection-reaction-diffusion problem, which we describe next.

3 Methods II: Numerical method To solve the fluid-poroelastic structure interaction and transport problems numerically, we use the diffuse interface method [10, 11]. Let χ denote the Heaviside function which equals one in Ω F and zero in Ω P . Let a phase-field function Φ: Ω → [0, 1] be a regularization of the Heaviside function such that Φ ≈ 1 in Ω F , Φ ≈ 0 in Ω P , and Φ smoothly transitions between these two values on a “diffuse” layer of width ϵ (see Fig 3). We suppose that dA ≈ |∇Φ|dV and . Using this notation, for functions F and f defined on Ω, we can write: where δ Γ is a Dirac distribution at the interface Γ. PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 3. Graphical representation of the diffuse interface approach in one dimension. Top: The phase-field function Φ. Bottom: The gradient of Φ used to approximate the location of the interface. https://doi.org/10.1371/journal.pcbi.1012079.g003 Using the approximations above, and , we can write the diffuse interface formulation of problem (8) as follows: Find such that for every , the following equality holds: (13) where we used the subscript to denote the “approximate” tangential component of a vector function at the diffused interface, defined, for a vector v, as follows: Tangent vectors that form a basis with ∇Φ can also be obtained directly from the phase-field function using the algorithm described in [12, 21]. We note that in order to write (13), the variables on each subdomain have to be extended onto the entire domain Ω. This procedure introduces singularities. For example, when the fluid velocity, which is defined on Ω F , is extended into the whole domain Ω, the corresponding integrand will be multiplied by Φ. Since Φ is zero in a large part of Ω P , the resulting linear system will have zero rows whenever this occurs, which will lead to a singular matrix. Therefore, we use the following regularization of Φ: (14) where β is a small positive number. Therefore, Φ ≥ β and 1 − Φ ≥ β. This regularization of Φ was used in the definition of problem (13). A phase-field method for the related Stokes-Darcy problems have been analyzed by the authors in [22] where existence of a weak solution was proved, and its convergence to the sharp interface solution was obtained. Because the concentration Eq (12) is already defined in Ω, we do not require regularization. In that case, β = 0. However, we use Φ to define the global velocity and diffusion: Numerical discretization Problems (13) and (12) are discretized in time using the Backward Euler method, and in space using a finite element method. Since the transport problem requires the fluid velocity, but the fluid-poroelastic structure interaction problem does not require any information from the transport problem, we first solve the fluid-poroelastic structure interaction problem until the final time T is reached. Then, we solve the transport problem. Let be the number of time steps for the fluid problem, T the final time, and Δt = T/N the time step. We define the discrete times tn = nΔt, for n = 0, …, N. We also denote by h a discretization parameter associated with the triangulation of Ω. For each h, we choose finite dimensional subspaces and over the triangulation . We use MINI elements [23] to approximate the fluid velocity and pressure, elements to approximate the Darcy pressure and the displacement, and elements to approximate the concentration. The fully discrete fluid-poroelastic structure interaction problem is given as follows: given for n = 0, …, N − 1, find such that for every , the following equality holds: (15) To solve the transport problem, we let be the number of time steps, and Δt t = T t /N t the time step. The fully discrete transport problem reads as follows: for n = 0, …, N t − 1, find such that for every , the following equality holds: (16) We note that according to (11), initially we have .

7 Discussion: Conclusions We devised a computational approach within the diffuse interface framework to explore the influence of scaffold architecture geometry on oxygen transport within biological scaffolds commonly employed in bioartificial organ engineering, with a specific focus on the bioartificial pancreas. To achieve this objective, we introduced a multi-physics model comprising a fluid flow component and an advection-reaction-diffusion model to analyze oxygen concentration within the scaffold. The fluid flow model incorporates the time-dependent Stokes equations coupled with the Biot equations, characterizing the behavior of a poroelastic medium representing the poroelastic hydrogel used in the scaffold design. We explored three biologically inspired scaffold architecture geometries: vertically drilled channels, branching channels, and a hexagonal geometry. From a computational standpoint, one of the primary challenges in addressing problems with multiple varying geometries lies in generating new computational geometries and designing appropriate matrices for spatial discretization. These matrices describe the unknown variables, such as fluid velocity, pressure, and oxygen concentration, across different computational domains. To streamline mesh generation in complex geometries, we introduced a diffuse interface approach in this study. In the diffuse interface approach, the unknown variables are defined across the entire scaffold domain, with the specific geometry of the channel network captured by redefining only the phase-field function. This simplification proves crucial not only for our current work but also for future endeavors, where we aim to develop a geometric optimization solver. This solver will simplify the generation of numerous channel geometries to optimize scaffold architecture. It is important to notice that a drawback of the diffuse interface method is the large size of the discretization matrices since the number of unknowns at the discrete level is doubled. Furthermore, to obtain higher accuracy, the mesh is commonly refined around the interface, which leads to a higher number of degrees of freedom. However, in simulations where the permeability is small, a fine mesh around the interface is also needed when a sharp interface method is used to accurately resolve the pressure gradient between the Stokes and Biot regions. Such mesh must align with the interface if a sharp interface method is used, while this is not required for the diffuse interface method. We demonstrated that the hexagonal geometry significantly outperforms both the branching channels’ network and the classical vertical channel geometries. Our analysis indicates that the superior performance of the hexagonal geometry stems from the relatively large angle between the dominant channel flow direction and the channel-hydrogel interface. This configuration results in a larger Darcy velocity, thereby facilitating enhanced advection-mediated oxygen supply to the transplanted cells. This study is significant because recent developments in hydrogel fabrication make it now possible to control hydrogel rheology [2, 5], utilizing the computational results to generate optimized scaffold architectures. Our future work includes the design of a geometric optimization algorithm for optimal scaffold architecture design. In case of geometric optimization, or just a moving domain problem, we expect it to be necessary to adapt the mesh when the phase function changes significantly. This can be done at each step, or every few steps, depending on how fast the phase field function evolves, by using the gradient of Φ.

Acknowledgments We would like to acknowledge the support by Dr. Shuvo Roy and his Bioartificial Pancreas Design Lab at UC San Francisco.

[END]
---
[1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1012079

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/

You are viewing proxied material from magical.fish. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.