Bone remodeling is a complex and critical process for maintaining bone health and preventing age-related diseases like osteoporosis. In silico models are powerful tools for understanding the cellular and protein properties that regulate bone remodeling, especially in aging conditions. For this purpose, we recently established a 3D multiscale micro-MPA model of trabecular bone remodeling using longitudinal in vivo data from PolgA mice, a mouse model of premature aging. Our model includes receptor-ligand kinetics, mechanomics, diffusion, and decay of cytokines that regulate cell behavior. We simulated trabecular bone remodeling in the sixth caudal vertebra of five mice over four weeks and evaluated the static and dynamic morphometry of the trabecular bone microarchitecture. We identified a configuration of the model parameters to simulate homeostatic trabecular bone remodeling, here named basal. Additionally, we produced anabolic, anti-anabolic, catabolic, and anti-catabolic responses with an increase or decrease in the levels of osteoprotegerin (OPG), receptor activator of nuclear factor kB ligand (RANKL), and sclerostin (Scl) produced by the osteocytes. We found that changes in OPG and RANKL levels were positively and negatively correlated with the average trabecular bone volume fraction (BV/TV) values, respectively. Conversely, changes in Scl levels produced small fluctuations in BV/TV in comparison to the basal state. Scl was deemed to be the main driver of equilibrium, while RANKL and OPG were shown to be involved in changes in bone volume fraction with potential relevance for age-related bone features. Moreover, bone remodeling is regulated by the interaction between different cells and tissues across many spatial and temporal scales. However, quantifying cellular and protein properties in vivo, particularly in aging conditions, presents challenges. To address this, we adapted the previously developed micro-MPA model to investigate the effects of aging and mechanical loading on bone remodeling in aging mice using in vivo micro-CT images and micro-FE models. Our results showed that the mechanoregulation of remodeling was altered due to aging, and the anabolic effect of an increase in bone mass with mechanical loading was no longer present in aged mice. We also demonstrated the importance of single-cell mechanomics in simulating bone anabolic response and the disruption of this response with aging. Higher numbers of osteoblasts and mineralization rate along with single-cell mechanomics successfully correctly simulated the anabolic response in young mice, while their disruption led to no anabolic response with loading in aged mice. In summary, our study shows the potential of micro-MPA modeling to advance our understanding of bone remodeling mechanisms in different physiological and age-related conditions and guide the development of new therapies for bone diseases.

Bone maintains its structure by repairing microdamage and functionally adapts to the mechanical environment by remodeling, where bone-resorbing osteoclasts and bone-forming osteoblasts play important roles through their interactions as well as mechano-biochemical couplings. An imbalance in the activities of osteoclasts and osteoblasts causes bone diseases such as osteoporosis, thereby increasing the risk of fracture. In these remodeling activities, osteocytes, which differentiate from osteoblasts and are embedded in the mineralized bone matrix, have been suggested to sense the mechanical environment and regulate the remodeling activities by communicating through cellular networks. Several signaling molecules and their transduction pathways in bone cells have been identified to understand the mechanism of bone metabolism and diseases. However, it is challenging to predict bone remodeling as a system because of the complexities of their interactions, mechano-biochemical couplings, changes in the mechanical environment, and the effects of drugs.

In this study, we focused on the role of osteocyte mechanosensing and communication and developed a computational simulation platform that predicts complex bone adaptation and microdamage repair by remodeling in silico based on the mathematical models of molecular and cellular behaviors. We believe that such in silico experiments of bone remodeling, combined with in vivo and in vitro experiments will be a beneficial tool to improve bone remodeling research.

First, the basics of these mathematical models and computer simulations will be summarized, and the in silico experiments will be discussed. Next, a mouse femur model will be used to demonstrate the adaptation of the trabecular bone by remodeling in response to applied forces, and the usefulness of the developed platform will be confirmed by reproducing the osteoporosis associated with reduced mechanical force and abnormal signaling conditions. Furthermore, the administration of several osteoporosis drugs that inhibit bone resorption and promote bone formation will be presented, showing the spatiotemporal behaviors of molecules and cells that are difficult to observe in vivo.

In addition to these mechano-biochemical couplings, microdamage accumulation and repair by remodeling are also involved to simulate the target remodeling. We have developed a mathematical model of bone microdamage according to continuum damage mechanics, which assumed that the mechanical properties of bone tissue and cell behavior depend on this variable. The validity of the proposed model will be discussed by simulating the mechanical adaptation for a single trabecula and the cancellous bone. These in silico experiments on bone remodeling will clarify the relationships between bone remodeling and bone microdamage accumulation. This simulation has the potential to improve the prevention and treatment strategies of bone diseases from biological and biomechanical viewpoints.

Can biology-inspired complexity be obtained without biochemical components? Can we replicate ubiquitous biological processes using only model physical building blocks like DNA-coated colloids that have simple but programmable interactions? The last decades have seen tremendous progress in understanding the self-assembly mechanisms that enable the formation of complex, sub-micron scale structures, but embedding these structures with bio-inspired functional behaviors remains a considerable challenge. Here, we demonstrate a scheme for transferring energy between two colloidal clusters, in analogy to ATP hydrolysis. By coupling the two clusters, we show how the one acting as a machine catalyzes a structural transition in the one acting as a fuel source, releasing energy that drives the machine into a higher energy structural state. The coupled system shows a significantly reduced mean-first passage time.
This work demonstrates that energy consumption, a fundamental and enabling biological process, can be replicated without complex biochemical reactions. In contrast, theories of active matter often focus on the effect of energy consumption, not on the mechanism itself. However, the mechanism is intimately connected to the type of physical phenomena that can result. In a next step, we extend the scheme by a third structural state in the machine. This allows us to convert energy into work by driving a net flux in the machine, which is not possible in equilibrium and requires a fuel source.

How thousands of microtubules and molecular motors self-organize in a dynamic fashion in spindles remains a very complex and largely unexplained phenomenon. However, recent data from large-scale electron tomography (EM) has enabled quantification of the positions, lengths and configurations of individual microtubules (MTs) in metaphase spindle in HeLa cells. We leverage the static ultrastructural data from the EM to quantitatively extract the density and orientational fluctuations of the individual microtubules around a well-defined coarse grained steady state. We then compare the equal-time correlations computed from the data with those derived from a coarse-grained active liquid crystal theory. Our preliminary analysis reveal that director-director correlations decay as ∼q^-2 consistent with the theory in which MTs are locally aligned through passive cross-linkers, molecular motors, and steric effects. The density-density correlations display ∼q^-4 scaling for large q indicating that MT sliding interactions alongside with diffusive-like motions of microtubules dominate at short wavelengths. In addition, from the EM data we can quantitatively extract the transverse fluctuation of an individual filament about its average tangent. The static structure factor associated with the transverse fluctuations decay as ∼q^-2 and their amplitude as ∼q^-4 revealing a thermal-like Fourier bending spectrum. However, the apparent persistance lengths inferred from the data is orders of magnitudes smaller than expected from thermal fluctuations alone, hinting at large nonthermal forces bending the microtubules. In summary, our work highlights how microtubule resolution data from EM in combination with simple coarse-grained theories can be used to describe mechanical behavior of spindles with measurable and interpretable parameters.

The piconewton forces generated by molecular motors carrying cargo along microtubules, or by microtubules polymerizing against the cell cortex or artificial boundaries, are sufficient to deform long microtubules. When microtubules are sparse in the cytoplasm, their deformations are disordered, characterized by high-frequency buckling and inducing only localized cytoplasmic flows. When the microtubules are instead arranged in a dense forest, the nature of the microtubule deformations and induced cytoplasmic flows can change dramatically, giving rise to long-range order and coherent flows. Using a combination of experiments, large-scale simulations of microtubules interacting hydrodynamically through a viscous fluid, and a coarse-grained theory for dense beds of filaments, we elucidate the mechanisms that underlie the self-organization of microtubule ensembles and their subsequent generation of cell-spanning rotation in two examples: cytoplasmic streaming in the Drosophila Melanogaster oocyte, and spontaneous rotation of artificially confined asters in Xenopus Laevis extract.

A computational approach will be presented for the estimation of the in vivo magnitude and spatial distribution of mechanical material properties of organs and other soft tissues from standard clinical imaging data. To this end, a new shape-based objective function, quantifying the difference between the measured and predicted tissue mechanical response is introduced. By utilizing shape, rather than deformation or related quantity, standard clinical imaging data can be utilized (e.g., without tagging) without further manipulation to approximate such measures from the images first. This objective function can then be implemented into an optimization-based inverse solution approach to estimate mechanical properties of tissues from initial and final shapes derived from clinical imaging data. This approach is an extension of prior work by the authors that used a standard discretized version of the Hausdorff distance as an objective function in an iterative approach to material parameter estimation [1]. A key component of the new inverse approach is constructing the geometry of the region of interest using a signed distance function. As such, a novel level-set framework is introduced for the objective function that is easily differentiable, and thus, able to be implemented into an optimization framework to estimate the material parameters that minimize the objective function with respect to a target shape with relative computational efficiency. A set of simulated inverse problems was used to evaluate the inverse solution estimation procedure based on estimating the passive elasticity of human ventricular walls from standard cardiac imaging data and corresponding hemodynamic measurements. In evaluating the results, emphasis will be placed on not just the accuracy of the material parameter estimates, but also on the computational expense (e.g., number of forward finite element analyses) required to approximate the target response. Various levels of heterogeneity will be considered in terms of the effect on solution accuracy and/or need for regularization. Additionally, sensitivity to model error will be explored.

REFERENCES

[1] Xu, J., Wong, T.C., Simon, M.A., and Brigham, J.C. A Clinically Applicable Strategy to Estimate the In Vivo Distribution of Mechanical Material Properties of the Right Ventricular Wall. International Journal for Numerical Methods in Biomedical Engineering, (2021)

For physics-based simulations, forward (UQ) and inverse uncertainty propagation still generate challenges for computationally demanding, real-world applications, like the ones appearing in many scenarios in Bioengineering. Two main obstacles pertain to a high stochastic dimension and many necessary forward solver calls caused by the statistical nature of the underlying algorithms.

A high stochastic dimension usually precludes the reliable use of surrogate-based approaches to reduce computational costs. Instead, we propose a Bayesian multi-fidelity procedure that exploits conditional distributions between two or more model fidelities in a small data regime (100 to 300 high-fidelity evaluations are necessary to train the conditional). The required online sampling is entirely shifted to the low-fidelity models. The latter can, for example, use simplified physics and coarser numerical discretization and only need to share a (nonlinear) statistical relationship with the high-fidelity model, giving a high degree of flexibility in selecting and creating such low-fidelity models. In both cases, the forward UQ and the Bayesian inverse problem, our approach results in an accurate posterior distribution, despite the inaccurate and noisy information the low-fidelity models provide. In the case of forward UQ, we call our approach Bayesian Multi Fidelity Monte Carlo (BMFMC), and for the inverse problem Bayesian Multi-Fidelity Inverse Analysis (BMFIA).

Especially for the Bayesian inverse problem, BMFIA brings further advantages. It shifts the necessary derivative evaluations to the low-fidelity model, allowing the analyst to exploit adjoint implementations for simplified physical models. We then perform Bayesian inference with efficient stochastic variational methods, which require solely evaluations of the lower-fidelity model. If the low-fidelity models cannot provide model gradients, we propose new developments in sequential black box variational inference schemes. The latter poses the variational optimization problem without needing model gradients. The formulation is more efficient than gradient-free sampling schemes, such as sequential Monte Carlo or Markov Chain Monte Carlo methods. It gives reliable posterior approximations for moderate stochastic dimensions (up to 50) where surrogate approaches are already prohibitive.

Another possibility for some challenging problems in Bioengineering is to drastically reduce the number of forward model calls via approximating the log-likelihood by a surrogate model. The advantage of such a log-likelihood approximation over the outputs of the forward model is that instead of a potentially highly dimensional model output, only a scalar value has to be approximated. To allow the scalability of the approach to higher stochastic dimensions, the training samples of the surrogate are adaptively selected, and the uncertainties of the surrogate are taken into account.

In this presentation, we will show the proposed methods' essential aspects and their application to different challenging problems in Biomechanics.

Changes in our body that occur due to illness and disease also bring mechanical changes, which are often only observed by clinicians poking and prodding tissues. These present opportunities for new diagnosis and monitoring technologies. In our work we have taken a tissue biomechanics approach to study the wound healing process. This area is particularly important due to the large costs that wounds place on healthcare resources (>£4.7 billion in the UK annually) and pain experienced by patients.

To assess how wound healing relates to mechanical changes, we used a mouse model of acute wound healing. We made 4 mm diameter wounds in the skin using a biopsy punch and then the skin was allowed to recover with mechanical and histological assessment at days 1, 3, 7 and 14 post-wounding. These correlated with the stages of wound healing – haemostasis, inflammation, proliferation, and remodelling. We excised the tissue and undertook tensile testing for global properties, image-based local strain assessment using Digital Image Correlation (DIC), and Optical Coherence Topography Elastography (OCE) to characterise the material changes in the skin. We correlated these to histology with H&E and Picrosirius red staining for wound-staging and collagen fibre characterisation respectively. Concurrently we developed a finite element model to aid the mechanistic interpretation of the data.

Our tensile testing results showed no discernible change in the hyperelastic moduli/coefficients when a wound was present in skin, compared to intact skin. This contrasted with literature which had shown the opposite. We believe that this is due to the more physiologically relevant strains we used to test the tissue, rather than the “test to failure” approach of others. To understand this further, we measured the local in-plane strains in the skin during tension. We observed re-organisation of the tissue strains with a compliant ring which reduces stress on the wound. This mechanism during healing appears to ensure that wounds are protected whilst healing. When skin was simulated by finite element models it became clear that a bulk hyperelastic approach did not sufficiently account for these changes, and fibrous models would be required. Our histology data showed how the fibrous components in the skin vary during healing, which indicates models will require time-variant structural changes. Our analysis of the wound mechanics by OCE helps identify sub-surface mechanical changes and we will share our progress on this.

The changes that occur in tissues during physiological changes, presents a substantial opportunity for our bioengineering community but we need to have both experiments and models that reflect the reality of tissue changes. Our work has shown that the experimentally derived data is most useful for computational model development only if both global and local structural changes are considered. Furthermore, the need for a model that adapts to the time-variant changes during a disease’s progression become more central. These, set against a backdrop of a varied population, increase the challenge of accurate physiological modelling but present a huge opportunity for the computational and experimental communities to work together.

Considerable research efforts have been devoted to the development of microstructurally-based anisotropic continuum constitutive models of biological soft tissues. These formulations typically rely on the definition of one or more vector fields representing the local orientation of biological fibres (i.e. collagen fibres) within a tissue. Methods to incorporate such a structural information into anatomically realistic micromechanical finite element (FE) models of soft tissues such as skin are still lagging behind, particularly when it comes to models aiming to capture the complex local three-dimensional (3D) architecture of the collagen fibres network.

In order to improve the predictive power of the next generation of biophysical models of soft tissues it is essential to develop robust methods and methodologies to seamlessly integrate the microstructural characteristics of the collagen network. Here, we developed such an approach combining high-resolution imaging of human dermis, fibre orientation image analysis, voxel-based mesh generation and micromechanical FE analyses.

Fresh full-thickness abdominal skin extracted during a cosmetic surgery procedure from a 39 years-old female Caucasian patient was commercially obtained (TCS Cellworks, Buckingham, UK). The skin sample was cut into 1 mm slices and processed for serial block-face scanning electron microscopy using a high contrast fixation protocol (SBEM Protocol v7_01_10, https://ncmir.ucsd.edu/sbem-protocol) and embedded in Spurr resin (Agar Scientific, Stansted, UK). The tissue block was then trimmed, glued onto an aluminium pin, sputter coated in gold/palladium and imaged in the serial block-face imaging system 3View® (Gatan, Inc., Pleasanton, CA, USA) mounted inside a Zeiss Sigma VP field emission scanning electron microscope with variable pressure mode (Carl Zeiss Microscopy GmbH, Jena, Germany) and imaged at 2.5kV. The acquisition was done at a sampling XY resolutions of 2500 x 2500 pixels (i.e. 8 nm pixel size) every 50 nm, resulting in 376 images. For development purpose, and due to the large size of the data set, the original image stack was cropped to generate a sub-stack with a 600×600×50 voxel dimension. The stack of 50 images was processed using a 3D orientation analysis algorithm based on the calculation of local structure tensors. A series of Python and Fortran programmes were written to generate a voxel-based hexahedral FE mesh with the option to assign fibre vectors at node, element or integration point level, import it into the FE package Abaqus® (Simulia, Dassault Systèmes, Johnston, RI, USA), code constitutive equations for invariant-based anisotropic hyperelasticity via Abaqus® UMAT subroutines, assign material properties to matrix and fibre phases, and run numerical analyses to study the micromechanics of the dermis. The main objective of our analyses was to quantify and understand the effects of faithfully capturing collagen fibre orientation on the homogenised micromechanics of collagen assemblies, and compare our technique to approaches considering a uniform orientation of collagen fibres, with or without fibre dispersion. Uniaxial/biaxial extensions and shear tests were simulated. The results of our numerical analyses did not only demonstrate the criticality of accounting for local fibre orientation but also the importance of distinguishing between the matrix and fibre phases as spatially independent domains within a tissue block.

Finite element models have been used for the last 50 years to investigate relationships between bone geometry, material properties, mechanical loading, and bone mechanoadaptation. FE models provide complementary information to experimental studies and clinical observations and can be used to explain experiments or design new experiments. FE models are also a critical part of the prosthetic design process to understand bone remodeling around the prosthesis. This historical perspective will cover the highlights of what we have learned from 50 years of bone modeling.

Bone is a complex hierarchically arranged structure which is constantly adapting to mechanical demands. To understand and model this so-called remodelling process, mainly three scales must be considered: the macroscopic (organ) scale, the mesoscopic (tissue) scale as well as the cellular scale. In trabecular bone, active osteoclasts and active osteoblasts remodel individual trabeculae at the cellular scale. On the so-called mesoscopic (tissue) scale, the size and shape of individual trabeculae changes, which on the macroscopic (organ) scale can be observed as a change in bone density. Macroscopic mechanical loading affects the stress/strain state in the mesostructure and thus influences the bone cell activities and the remodelling process.

To capture these multiscale interactions, the aim of this work is to develop a three-scale approach to bone remodelling including mechanical feedback that combines the advantages of already established models to describe remodelling at the macroscopic, the mesoscopic and the cellular scale. Therefore, a continuum approach based on the theory of open system thermodynamics is used to represent the remodelling process on the macroscale. Approaches on the macroscopic scale are usually based on average material properties and locally adjust bone density to mechanical use following Wolff's law. Our macroscopic approach can also account for non-linearity, is efficient and numerically stable. Moreover, it is easily extendable to other aspects such as anisotropy, age dependence and biological availability of nutrients and hormones. However, it does not explicitly consider the very irregular trabecular structure of cancellous bone. On the mesoscale of our approach, the trabecular structure is modelled in a simplified way as an ideal truss network, where each trabecular is represented by a truss element. This approach requires much less computational effort than pixel- or voxel-based models, which are able to represent real geometries but are considered inefficient in a multi-scale context. On the microscale, a bone cell population model is used to capture the coupled activities of bone cells in the remodelling process. The main focus is first on the development and verification of a remodelling algorithm that captures all three scales. The remodelling algorithm is tested on benchmark problems of mechanical overload and disuse and applied to more complex geometries.

Mechanical loading is essential for maintaining homeostasis and driving adaptation of the skeleton. Disruption of this mechanically-driven bone remodeling process has been linked to bone loss and increased fracture risk. The etiology of bone pathology is complex, particularly in cases involving metabolic bone diseases, and no current serological biomarkers directly reflect bone mechanical quality. Rather, indirect measures like bone mineral density, clinical risk prediction tools, or visual interpretation of imaging are used to drive treatment decisions. As a result, the prognostic limitations of most clinical assessments reinforce patterns of care that tend to be conservative and reactive waiting for clinical pathology to present unambiguously – rather than proactive and preventive. Despite these challenges, innovation in medical imaging and increasingly accessible computational power have made it possible to computationally measure the mechanical properties of bone via image-based finite element analysis. Moreover, longitudinal clinical studies have also shown that the combined use of high-resolution peripheral quantitative computed tomography (HR-pQCT) and micro-finite element analysis can be used to measure load-driven bone remodelling in vivo. This powerful combination enables the assessment of bone formation and resorption (i.e., remodeling) and quantification of the relationship between this remodeling and mechanical loading (i.e., mechanoregulation) at the microstructural level over time.

Despite the promise of such advanced imaging and computational methods in clinical studies, there remain technological challenges to widespread implementation. These challenges include variations in image quality due to image noise and patient motion, which confound the precision and reproducibility of bone density, geometry, mechanical, cortical, and trabecular structure measurements. Additionally, segmentation protocols and image processing pipelines often require manual input or augmentation, prohibiting the high-throughput analysis that would be required for clinical implementation. Finally, to push such imaging and computational methods for mechanobiological bone remodeling studies from bench (supercomputer) to bedside (clinical computing) in the hospital environment these methods need to be repeatable and validated.

With this in mind, this talk will (1) explore the challenges associated with the assessment of dynamic bone morphometry in vivo; (2) link bone microarchitecture and functional outcomes using the finite-element method in patient populations with metabolic bone disease such as osteoporosis and diabetes mellitus; and (3) evaluate the potential for such in silico modelling to serve as a clinical tool for monitoring changes in bone health over time. With continued evolution of technology and best practices, the utility of bone mechanical biomarkers using image-based finite element analysis will undoubtedly increase, revolutionizing standard of care in bone health.

Understanding how the loss in the mechanical environment is associated with bone loss is challenging. Endocortical bone loss due to the disuse of muscles has been recently reported in the literature[1], which shows that bone loss was focused on the posterior surface at the distal section; however, it shifted to the anterior lateral surface at the proximal end. Interestingly, it conveys that bone loss is site specific along the bone length. However, the cause of bone loss in disuse conditions is unknown. In this work, we aim to decipher such an underlying mechanism.

Recently load-induced fluid flow inside the lacune canaliculi network (LCN) has caught the attention of researchers as a primary stimulus. In loaded bone, fluid flow will be higher at the endocortical surface and negligible at the periosteal surface due to its permeable and impermeable nature, respectively. It seems reasonable to assume that in the case of a disuse condition, the maximum loss of fluid flow will occur at the endocortical surface.

Hence, the present work explores the loss of fluid flow related to disuse and aims to develop a mathematical model that can correctly predict the site of bone loss. We use the loss of dissipation energy density due to disuse as a stimulus as it follows the trend of loss of fluid velocity. Accordingly, we hypothesized that the site of maximum loss of dissipation energy density due to fluid flow corresponds to the site of maximum loss of bone tissue. To test our hypothesis, a poroelasticity-based finite element model of a simplified geometry was subjected to the physiological loading condition. The fluid velocity computed from the above analysis is used to calculate the dissipation energy density (φ_phys). Simultaneously, we assume that the disuse of bone causes no fluid flow, and the corresponding dissipation energy density is denoted as (φ_disuse). Therefore, the dissipation loss can be written as φ_phys- φ_disuse: this dissipation loss and experimental bone loss data reported by Ausk et al.[1] are coupled to develop a novel mathematical model that predicts the site of cortical bone loss with spatial accuracy.

The developed mathematical model first predicts bone loss given in the literature at distal (p value= 0.89, Watson U² test) and midshaft (p value= 0.78, Watson U² test). However, it overestimates bone loss at the proximal section. Based on the availability of experimental data, the developed model can be extended to predict other disuse conditions such as spaceflight, prolonged bed rest, etc.

[1] B. J. Ausk, P. Huber, S. L. Poliachik, S. D. Bain, S. Srinivasan, and T. S. Gross, “Cortical bone resorption following muscle paralysis is spatially heterogeneous,” Bone, vol. 50, no. 1, pp. 14–22, Jan. 2012, doi: 10.1016/j.bone.2011.08.028.

Sepsis, a life-threatening disease caused by an unregulated immune response to infection, is a significant global health concern, contributing to approximately one fifth of all worldwide deaths. The chances of a patient surviving sepsis decrease by 8% for every hour before treatment is started, making rapid diagnosis critical to improving patient outcomes. Current sepsis diagnosis methods may take 24 to 48 hours, highlighting the need for faster alternatives. To address this challenge, a novel sepsis detection technique (Guillou et al. 2021) has been developed that evaluates the behaviour of leukocytes as they traverse a microfluidic cross-slot junction. The cell behaviour, specifically the degree of cell stretching due to elongational stresses in the junction, and the rate at which the cell oscillates when interacting with a vortex present in the junction outlets, are measured and used in statistical models to interpret the results. This approach enables the rapid and accurate determination of an individual's risk of developing sepsis within 10 minutes, significantly improving patient outcomes.

Despite the technique's proven ability to ability to correlate leukocyte behaviour in the microfluidic junction with patient sepsis risk, there is no theoretical understanding that connects the measured cell behaviour with the pathophysiological changes that differentiate septic and healthy cells. Previous studies have simulated the behaviour of rigid spheres in the microfluidic cross-slot and have shown that the device is highly sensitive to both the initial position and the size of particles entering the junction. (Kechagidis et al. 2022) However, leukocytes have a complex mechanical structure which influences the way in which they interact with fluids, especially when deforming due to fluid stresses. This means that to build a full picture of the behaviour of cells within the device, considerations of the leukocyte mechanical structure must be included.We have developed a detailed leukocyte model capable of simulating cell membrane viscoelastic and non-linear behaviour, cytoplasm viscosity, and non-homogeneous considerations such as cell nuclei. The cell model employs the finite element method for tracking cell membranes, coupled via the immersed boundary method to an extensively validated custom lattice-Boltzmann fluid solver.

In this talk, I will present the results of a series of simulations that examine the sensitivity of sepsis detection metrics in the microfluidic device to the parameters employed in the leukocyte model, such as the viscosity or stiffness of cells and their nuclei. By analysing the interaction of the flow field with the leukocyte and its internal structure, I will discuss the mechanisms underlying this sensitivity.

References
Guillou, L., Sheybani, R., Jensen, A.E., Di Carlo, D., Caffery, T.S., Thomas, C.B., Shah, A.M., Tse, H.T.K., O’Neal, H.R., 2021. Development and validation of a cellular host response test as an early diagnostic for sepsis. PLoS ONE. https://doi.org/10.1371/journal.pone.0246980

Kechagidis K., Owen B., Guillou L., Tse H., Di Carlo D., Krüger T., 2022. Numerical investigation of the dynamics of a rigid spherical particle in a vortical cross-slot flow at moderate inertia. BioRxiv https://doi.org/10.1101/2022.12.19.520995

Elucidating the lymph fluid pumping mechanism is of crucial significance towards a better understanding of the lymphatic system related diseases, such as lymphedema and cancer, and subsequently the development of more efficient treatments. To this end, we develop numerical framework to model the interplay between the lymph fluid flow, the contraction-relaxation of the lymphatic vessel walls and the lymphatic two-leaflet valves opening-closing. Our numerical method accounts for the biochemical signals of calcium ions (Ca²⁺) and the nitric oxide (NO) regulating the vessel contraction-relaxation. We use the lattice Boltzmann method to compute the lymph fluid flow and the chemical species mass transport, a spring network model for the lymphatic vessel walls and valves, and the immersed boundary method to implement the two-way coupling of the fluid-structure interaction. Results on the effects of various fluid, geometrical and mechanical properties on the mechanical performance of the valves, the vessel and the lymphatic pumping mechanism will be presented.

Thrombotic deposition is a major consideration in the development of implantable cardiovascular devices. The process of thrombosis is governed by both biochemical and mechanical considerations. The mechanical contribution has traditionally been thought to involve platelet activation due to sufficiently elevated shear stress levels in blood. Recent experimental evidence, however, suggests that beyond the effect of shear stress alone, localized changes in the blood shear rate, i.e. shear gradients, play a critical role in thrombogenesis by controlling the process of platelet deposition. The goal of the present work is to develop a predictive computational model of platelet plug formation that incorporates for the first time the effects of shear gradients and that can be used to assess the thrombotic burden of cardiovascular devices. To this end, we have developed a comprehensive model of platelet-mediated thrombogenesis which includes platelet transport in blood flow, platelet activation and aggregation induced by both biochemical and mechanical factors, and the kinetics and mechanics of platelet adhesion. Moreover, we also consider the effect of thrombus growth on the local fluid dynamics and how these alterations in the fluid dynamic environment feed back into the process of thrombogenesis.

The two-dimensional computational model was developed using the commercial multi-physics finite element solver COMSOL 5.6. The biochemical component of the model can be described by a set of coupled convection-diffusion-reaction equations that account for resting platelets, activated platelets, and various pro- and anti-thrombotic chemical species. Platelet adhesion to the blood vessel surface was modeled via a flux boundary condition. By incorporating a moving mesh for the blood vessel wall into the model, thrombus growth and consequent alterations in blood flow were modeled. The mechanical component of the model includes the induction of platelet activation by sufficiently high levels of shear stress and the promotion of platelet deposition in zones of negative shear gradients. Two physiologically and pathologically relevant scenarios were studied. The first scenario involves stenoses of varying severity where the mechanical contribution includes platelet activation in the contraction zone and platelet deposition in the expansion zone downstream of the stenosis. The second scenario involves a bifurcation where local disturbances in the flow field dictate the localization and extent of the computed thrombi.

The results demonstrate the model’s ability to provide the spatial and temporal evolution of a platelet plug within a flow field. The computed platelet plug size evolution was validated against experimental data from the literature. The results confirm the importance of considering both the chemical and mechanical contributions to platelet aggregation and underscore the importance of accounting for the effects of shear gradients. The developed model represents a potentially useful tool for the optimization of the design of the cardiovascular device flow path.

Interendothelial slits in the spleen fulfill the major physiological function of continuously filtering red blood cells (RBCs) from the bloodstream to remove abnormal and aged cells. To date, the process of extreme deformation and passage of 8 µm RBCs through 0.3-µm wide slits remains enigmatic. Should the slits increase their caliber during RBC passage as sometimes proposed in the literature? The values of the mechanical quantities selected during the passage of RBCs in the splenic slits and the associated dynamics, such as deformation mechanisms and transit times, remain poorly known today. Recent numerical and/or theoretical approaches have suggested that the spleen may play an important role in defining the surface area-to-volume ratio of the RBCs circulating in the microvascular system, but, so far, numerical approaches are not quantitatively validated by experiments and no experimental direct observation of RBCs flowing in splenic-like slits supports this hypothesis.

Here, we couple a unique in-vitro microfluidic technique to a multiscale in-silico RBC model that enables a quantitative approach of the mechanisms of passage of RBCs through interendothelial slits. The in-vitro technique allows the observation of the dynamics of passage of RBCs in slits of submicron width under tunable external stresses. The in-silico RBC model is implemented in dynamic and quasi-static versions. The dynamic version is integrated with a boundary integral simulation of surrounding flows to resolve the full fluid-cell interactions during this passage process, while the quasi-static version is done in commercial software ABAQUS. Agreement between experiments are simulations is remarkable.

We show that RBCs are capable of amazing extreme deformations allowing them to pass through rigid slits as narrow as 0.28 μm under a pressure drop of 500 Pa at body temperature, but not at room temperature. To achieve this tour de force, they must meet two requirements. Geometrically, the surface area-to-volume ratio of individual cells must be sufficient to form two tether-connected equal spheres. Mechanically, they must be able to locally unfold their spectrin cytoskeleton inside the slits. In contrast, activation of the mechanosensitive PIEZO1 channel is not required. The RBC transit time through slits scales with in-slit pressure drop and slit width to the -1 and -3 power, respectively. This transit dynamics is similar to that of a Newtonian fluid in a 2D Poiseuille flow, thus showing that it is controlled by the RBC cytoplasmic viscosity. We quantitatively predict transit times as a function of the cell mechanical properties and external parameters: pressure drop, slit size, and temperature. Altogether, our results clearly show that filtration through submicron-wide slits is possible without further slit opening. Furthermore, our coupled experimental/simulation approach is quantitative and addresses the critical need for in-vitro evaluation of splenic clearance of diseased or engineered RBCs for transfusion and drug delivery.

A.Moreau, F. Yaya, H. Lu, A. Surendranath, A. Charrier, B. Dehapiot, E. Helfer, A. Viallat, Z. Peng, bioRxiv 2023.01.10.523245; doi: https://doi.org/10.1101/2023.01.10.523245

Model simulations indicate that the response of growing cell populations on mechanical stress follows the same functional relationship and is predictable over different cell lines and growth conditions. We developed an Agent-Based hybrid model strategy in which cells are represented by coarse-grained individual units calibrated with a high resolution cell model and parameterized by measurable biophysical and cell-biological parameters. After parameter calibration from experiments with mouse colon carcinoma cells growing against the resistance of an elastic alginate capsule, the model adequately predicts the growth curve in i) soft and rigid capsules, ii) in different experimental conditions where the mechanical stress is generated by osmosis via a high molecular weight dextran solution, and iii) for other cell types with different growth kinetics from the growth kinetics in absence of external stress. Our model simulation results suggest a generic growth response of cell populations upon externally applied mechanical stress, as it can be quantitatively predicted using the same growth progression function.

During early embryonic development, cells organize into multicellular arrangements of epithelial tissue that continuously grows and deforms without losing tissue integrity. The embryo shape changes depend on the interplay between stress acting on or within the epithelial sheet and the ability of the epithelium to remodel, i.e. tissue rheological properties. Cell divisions appear to play a fundamental role in this process, since it not only drives tissue growth, but can also affect tissue mechanical stress and tissue remodeling by cell mitotic rounding [1], and the forces generated during cytokinesis. In turn, mechanical forces regulate the rate and orientation of cell division through mechanosensitive feedback processes [2]. A well-known mechanism is contact inhibition of proliferation where high cell density surpresses cell division. Defects in these feedback loops can result in uncontrolled tissue growth and loss of tissue integrity and epithelial structure, commonly observed during tumor growth.

Understanding the dynamics and mechanics of epithelial tissues during embryonic development requires computational modeling from the single-cell perspective. Therefore, we developed a 3D deformable cell model that accurately describes cell mechanical properties, cell-cell interactions and its accompanied shape changes before and during cell divisions. In this model, a cell is represented as a fluid-filled, viscous vesicle under active tension with adhesive cell-cell interactions analogous to a bubble inside a foam. Finally, cytokinesis is driven by an elastic contractile furrow ring [3].

In this work, we implement a mechanosensitive feedback loop between cell growth and division, and mechanical stress perceived by the individual cell inside the tissue. Next, this model is applied to study the tissue steady-state (homeostatic) mechanical state such as tissue stress, cell density,… in function of single-cell mechanical properties such as cell-cell adhesive interactions, cell cortex contractility and cell growth pressure. Doing so, allows us to unravel the interplay between cell growth, division and epithelial integrity.

References

[1] Petridou N., et al. (2019). Fluidization-mediated tissue spreading by mitotic cell rounding and non-canonical Wnt signalling. Nature Cell Biology. Vol 21

[2] Godard, B., Heisenberg C. (2019). Cell Division and tissue mechanics. Current Opinion in Cell Biology

[3] Cuvelier, M., et al. (2022). Stability of asymmetric cell division under confinement: A deformable cell model of cytokinesis applied to C. elegans development. bioRxiv

One of the main challenges in improving the efficacy of conventional chemotherapeutic drugs is that they do not reach the cancer cells at sufficiently high doses, while at the same time affecting healthy tissue and causing significant side effects and suffering in cancer patients. To overcome this deficiency, magnetic nanoparticle-based drugs have emerged as a promising approach to achieve more specific tumour targeting. Nanoparticles are nano-sized organic or inorganic materials that can be equipped with various biological and medical functions and designed with different physicochemical properties: a chemotherapeutic agent can be encapsulated in the nanoparticles or attached to their surface, and magnetic nanoparticles can be directed to the target tissue by applying an external magnetic field. However, the geometry and position of the magnet are critical to the effective delivery of the drug to the tumour.

We therefore present a computational model for the magnetic targeting of nanoparticles that are designed to deliver therapeutic agents to the tumour: our model includes the interstitial flow around the tumour, the magnetic forces that guide the nanoparticles, and the transport within the tumour. We show how a model for the transport of magnetic nanoparticles in an external magnetic field can be integrated with a multiphase tumour growth model based on porous media. Tumour spheroids as model systems are on the scale of a few hundred micrometres, whereas nanoparticles are several orders of magnitude smaller. While it is possible to computationally study the movement of individual particles, the computational burden is enormous when thousands of particles are involved at the scale of an entire tumour, and we hence use a homogenised approach at the macroscale.

Such a computational model allows studying the influence of the geometry and position of the magnet on the accumulation of nanoparticles in the tumour, which is essential for successful cancer therapy. The ultimate goal of any cancer therapy is to optimise the effectiveness of the treatment strategy, and to this end we also show how probabilistic Bayesian methods can be used to optimise magnetic drug targeting. Our entire approach is based on the underlying physical mechanisms, not just on big data or artificial intelligence. It can therefore provide crucial insights into mechanisms that cannot be studied conclusively in vivo or in vitro alone. Model development, in vivo or in vitro experiments, and clinical trials must complement each other in a feedback loop to achieve more personalised treatment strategies that improve therapeutic outcomes and limit adverse effects for cancer patients.

A short remaining life expectancy and high mortality characterise brain tumours as a particularly dangerous disease. Simulations of the relevant processes of tumour growth and regression in brain tissue are realised by embedding a continuum-mechanical model in the framework of the Theory of Porous Media (TPM) into a data-integrated workflow. This workflow is based on suitable patient-specific data, that is basically available, e. g. from magnetic resonance images (MRI). Preparing the data by a set of tools, for example using a convolutional neural network in a shape of an U-Net, result in the segmented position and composition of the tumour and provide the referential geometry of an initial boundary value problem (IBVP). Furthermore, relevant information, e. g. about heterogenities or flow properties, are collected by image-processing tools. A more cost effective surrogate model based on the ratio of the composition of the tumour compartments is developed and calibrated with simulations of the TPM model. These modularly arranged components of the developed data-integrated approach are processed using the finite-element framework FEniCS and allow to study relevant clinical questions.

Cancer cell mechanics dictates cellular force generation, intra- and inter-cellular force transmission, cell morphology, cell migration and cell metastatic potential. The mechanics is a function of organization and dynamics of several cytoskeletal elements such as actin, microtubule and intermediate filaments, molecular motors interacting with these filaments, filament cross-linkers and other associated proteins that guide filaments polymerization, branching and depolymerization. These elements are further influenced by the activity of various enzymes, metabolites and small molecule effectors. The biochemical signaling networks upstream of a cell’s cytoskeletal elements can be complex with competing interactions, vary based on cellular genotype and phenotype, and consequently drive diverse mechanical behaviors in cells in response to microenvironmental signals. We map this extensive cytoskeletal signaling network, develop a new Boolean-Hybrid-Modular model to solve this signaling network and use it to predict changes in cytoskeletal elements and their organization in response to external and internal signals. We further connect these changes in cytoskeletal elements with existing cell-mechanics based models of traction force generation and migration to understand how a confluence of various micro-environmental signals and cell phenotypes drive disparate cancer cell behaviors.

Bacteria predominantly reside in multicellular communities, called biofilms, rather than as independent planktonic organisms. These biofilms are encompassed by a self-produced extracellular matrix that provides mechanical support and serves as a barrier for external stressors such as antimicrobial treatments. Hence, biofilms are notoriously hard to treat and eradicate and consequently confer a large societal and medical toll. Although the overall protective effects of extracellular matrix production are well-known, it remains unclear how the production of matrix at single-cell level can shape the structure, and hence the protection, of the bacteria and matrix. Using a combination of individual-cell based modeling and in vitro single-cell level characterization of Salmonella biofilms, we aim to elucidate the emergence of biofilm structure in function of matrix productivity. In silico, we model growth-driven biofilm formation where matrix production is represented through the modulation of cell-cell adhesion and friction. We find that increases in cell-cell adhesion result in more aligned and densely packed biofilms, whereas increases in cell-cell friction decrease orientational order and density. As such, increases in matrix production can potentially result in more porous biofilms, and, contrary to the current paradigm, less protection against antimicrobial treatment.

Cardiovascular diseases are the main cause of death worldwide. The myocardial bride is one of the congenital abnormalities while part of the coronary artery is placed under the myocardium instead of resting on top of it. The myocardial bridge can be asymptomatic during adult life or can cause many cardiovascular consequences such as angina, myocardial ischemia, acute coronary syndrome, left ventricular dysfunction, arrhythmia, and even sudden cardiac death. Therefore, the investigation of the presented abnormality can significantly influence its understanding.

The presented research focuses on the numerical modelling of the blood flow through the myocardial bridge section. Computer modelling is a useful tool to assess the risk factors influenced by blood flow disturbances in arteriosclerosis disease development. The mentioned factors are mainly connected with the low shear stress affected blood vessels accompanied by high oscillations of the blood flow in specific regions. In the case of the myocardial bridge, the blood vessel movement and its contraction can significantly influence the results of the flow calculations. Therefore, numerical models should include the simulated arteries' geometry changes, which is a challenging model development. It can be realized by the dynamic mesh approach accounting for the dynamic discretization of the investigated domain during the calculations. The presented research describes two methods of mesh motion defining in the commercial CFD software (ANSYS Fluent) extended by additional user defined functions. The first is based on the shape reconstruction from the camera records collected in in-vitro conditions on the blood vessel phantom. Results from the model and experimental campaign are compared including shape, pressures and flow in the blood vessel phantom. The pressure range was in this case at the level of approx. 120/80 mmHg. The consistency between those results was at a satisfactory level. The second mesh motion deformation was based on the patient's medical images collected in in-vivo conditions and processed by auxiliary software (ANTS) performing advanced image transformation based on diffeomorphism. The obtained results from those two methods pointed out the validity of the dynamic domain definition. The relevant volume change of the vessel during the heart cycle also influences the assumed boundary conditions introduced to solve governing equations covering specifically the mass conservation principle. The dynamic shape reconstruction in the CFD model was compared with the instant segmentation data showing satisfactory agreement.

CFD is a noninvasive technique that could be applied to clinical practice to predict the blood flow in the coronary arteries. Including the vessel shape dynamics from medical images in the CFD models increase those models' accuracy.

Acknowledgements: This research is supported by National Science Centre (Poland) project No. 2017/27/B/ST8/01046 and project No. 2019/34/H/ST8/00624. This help is gratefully acknowledged.

Coil embolization is performed to treat intracranial aneurysms. During the follow-up after coil embolization, a recurrence of the aneurysm called recanalization may occur due to blood inflow. Therefore, it is expected that risk prediction of aneurysmal recanalization will enable effective treatment planning. Recently, computational fluid dynamics (CFD) has been conducted to identify the risk factors of recanalization. Although most of the previous studies assumed that blood is a Newtonian fluid, actual blood is a non-Newtonian fluid. Especially, the shear rate in coil embolized aneurysms becomes so low that the effect of non-Newtonian fluid is considered to be significant. In this study, we conducted CFD and statistical analysis to investigate the effect of non-Newtonian fluid properties on the prediction of aneurysmal recanalization. Our target of the analysis was basilar tip aneurysms. We identified 31 aneurysms including 4 recanalized aneurysms and 27 stable aneurysms. In this study, we defined recanalization as the case that recanalized less than 2 years after the first coil embolization. On the other hand, aneurysms that remained stable for more than 2 years after the embolization were defined as stable. We applied Newtonian and non-Newtonian blood models for these cases in the CFD analysis. We modeled the non-Newtonian blood using a modified Casson model based on the results of viscosity measurements on 12 patients with aneurysms. A porous media model was applied to model the embolized coil mass. We defined 33 hemodynamic and 3 morphological parameters. For each hemodynamic parameter, the value before (BF) and after (AF) coil embolization was obtained, and their reduction ratio (RR) was also calculated. We focused in particular on the V_Neck and the V_Dome, which indicate the dimensionless velocity on the neck plane in the aneurysm, respectively. In addition, the spatial mean, maximum, and minimum values were denoted as Ave, Max, and Min, respectively. Univariate logistic regression analysis was performed to compare each parameter between recanalized and stable cases. In addition, ROC analysis was also performed for parameters that showed statistically significant differences (P < 0.05). The sensitivity, specificity, and area under the curve (AUC) were evaluated among the parameters. As a result, by considering non-Newtonian fluid properties, V_{Neck_Ave_RR}, V_{Neck_Max_RR}, and V_{Dome_Ave_RR} newly indicated statistically significant differences between the recanalized and stable cases (i.e., the parameters that showed statistically significant differences in the Newtonian fluid also showed statistically significant differences in the non-Newtonian fluid). In addition, although the highest sensitivity (=0.75) and AUC (=0.84) were obtained with the V_{Dome_Ave_RR} when we assumed Newtonian fluid, the higher sensitivity (=1.00) and AUC (=0.89) were achieved when the non-Newtonian model was applied. The results showed that considering the non-Newtonian properties increased the sensitivity and AUC of recanalization prediction. In conclusion, non-Newtonian fluid properties changed the parameters that showed statistically significant differences between recanalization and stable cases. Furthermore, the V_{Dome_Ave_RR} obtained the highest accuracy in all researched parameters assuming Newtonian and Non-Newtonian fluid. The present result implies that the introduction of non-Newtonian fluid properties can improve the prediction accuracy of recanalization in coil embolized aneurysms.

Understanding cardiac sympathovagal balance and its relationship to the degree of reinnervation after heart transplantation has been a longstanding challenge in medicine. While autonomic markers derived from cardiac rhythms are often used to assess cardiac sympathovagal balance, their relationship to reinnervation has not been fully explored. In this study, we used a mathematical model of the human cardiovascular system and its autonomic control to explore the influence of varying levels of vagal and sympathetic cardiac reinnervation on autonomic cardiac markers.

To investigate this relationship, we applied a mathematical model of the human cardiovascular system and its autonomic control. The model integrates the chronotropic and inotropic effects of the arterial baroreflex and pulmonary stretch reflex, which are the main contributors to heart rate variability. The work focused on markers that are commonly used to assess cardiac sympathovagal balance, including resting heart rate, bradycardic and tachycardic heart rate response to the Valsalva maneuver, root mean square error of normalized RR-intervals (RMSDD), high-frequency (HF) power, low-frequency (LF) power, and total power of the heart rate variability spectrum. To evaluate the strength of the relationship between the level of cardiac reinnervation and the respective markers we calculated Spearman's rank correlation coefficients.

Results showed that for assessing vagal cardiac reinnervation levels above 20%, resting heart rate, RMSDD, and total spectral power may be equally suitable as the commonly used measure of HF power. The strength of the correlation between these markers and vagal reinnervation was found to be very high, with correlation coefficients of ρ=0.99, ρ=0.97, and ρ=0.89, respectively (all p<0.05). Concerning sympathetic reinnervation, simulations suggest that LF/HF-ratio and tachycardic response to the Valsalva maneuver may be more suitable than the regularly used measure of LF-power. The strength of the correlation between these markers and sympathetic reinnervation was high: ρ=0.88 and ρ=0.84, respectively (both p<0.05). These findings suggest that there are differences in the performance of cardiac autonomic markers in assessing vagal and sympathetic reinnervation and that some markers may be more suitable than others depending on the level of reinnervation.

The developed mathematical model can provide critical insights into the genesis of autonomic cardiac markers and their relationship to cardiac reinnervation. Furthermore, it suggests strategies for designing and interpreting future clinical studies, which would then provide clinical evidence for the current findings and allow more accurate evaluation of cardiac (re)innervation levels after heart transplantation. Overall, the findings of this model study might have important implications for the assessment and management of cardiovascular disease and highlight the potential of mathematical models to enhance our understanding of complex physiological systems.

This work was funded by the European Project H2020-EU.1.2.2. “A neuroprosthesis to restore the vagal-cardiac closed-loop connection after heart transplantation, NeuHeart” (Grant agreement ID: 824071).

Objectives: A standard treatment for severe Aortic Stenosis (AS) is surgical valve replacement with biological prostheses. Fluid-structure-interaction (FSI) simulations are a powerful tool for studying hemodynamic valve behavior and thus optimize function, geometry and durability particularly of prosthetic aortic valves. In this study, an FSI simulation of a bioprosthetic bovine pericardium aortic valve was established, which incorporated an anisotropic hyperelastic material model for the valve leaflets.

Methods: A valve-specific geometric model from a µCT scan (voxel size 31µm) of the 21 mm Avalus™ Bioprosthesis (Model 400, Medtronic, Ireland) was obtained. An anisotropic hyperelastic material model was used, a first order Ogden model for the isotropic matrix and an exponentially-based strain energy function for embedded hyperelastic fibers representing collagen fibers in circumferential leaflet direction. The valve structure was then coupled with a fluid domain consisting of a tube with 24mm diameter as inlet, the sinuses of Valsalva and a tube with 30mm diameter as outlet. Blood was assumed to be a Newtonian fluid. The inlet pressure boundary condition was chosen to represent a clinically determined transvalvular pressure gradient during a cardiac cycle, while the outlet pressure boundary condition was set to zero. The valve hemodynamic behavior was analyzed by determining the valve’s effective orifice area, leaflet fluttering dynamics during systole and fluid velocities.

Results: The composite valve model was successfully implemented in the FSI simulations performed on the Vienna Scientific Cluster supercomputer. The simulated effective orifice area at mid-systole was 1.6 cm². Systolic leaflet fluttering – which is associated with fatigue and failure of bioprosthetic leaflets over time – was most pronounced at beginning of systole, with amplitudes up to 2.5 mm. However, during the course of systole, the leaflets stabilized, resulting in smaller fluttering amplitudes but higher fluttering frequencies (up to the simulation output sampling rate of 200 Hz). Peak velocities in the free jet during systole went up to 3.0 m/s.

Conclusions: A hyperelastic anisotropic leaflet material was successfully integrated in an aortic valve model and used for FSI analysis. The model is being refined and validated with experimental flow tests. After validation the influence of differently sized valves on the valve hydrodynamics will be investigated together with correlations of simulation results with clinical outcomes.

An intracranial aneurysm is a cerebrovascular disease that swells of bulges abnormally in the wall of intracranial arteries. Rupture of aneurysms cause subarachnoid hemorrhage with high mortality. Unruptured aneurysms occasionally have Thin-Walled Regions (TWRs) where the wall thickness is thinner than the surrounding areas. Although TWRs in intracranial aneurysms have a risk of rupture, imaging modalities cannot evaluate the thickness of the aneurysm wall. The surgical treatment will be able to perform safely by identifying TWRs before the treatment. TWRs have been related to hemodynamic factors. Although Computational Fluid Dynamics (CFD) analysis has been applied to identify the hemodynamic characteristics around TWRs, few previous studies have focused on the difference in hemodynamic factors between TWRs and non-TWRs by quantitative definition of TWRs. The purpose of this study is to investigate the hemodynamic factors involved in TWRs by comparing the results in TWRs and non-TWRs by CFD and image analysis. We identified 100 aneurysms (middle cerebral artery: 78, anterior cerebral artery: 20, internal carotid artery: 2) treated with craniotomy and clipping. Unsteady CFD analysis was performed with the inlet boundary condition imposed on the mean mass flow pulsations measured from healthy adults, the outlet boundary condition fixed at a static pressure of 0 Pa, and the no-slip wall boundary. We evaluate the pressure difference (PD), wall shear stress (WSS), and wall shear stress divergence (WSSD) on the aneurysm wall. All parameters are normalized by the dynamic pressure at the aneurysm inlet. Since TWRs generally indicate intense red, the comprehensive Red (cR) value was defined by the RGB color model to evaluate the redness. The cR value was calculated for each pixel of the intraoperative images, and TWRs were defined using the cR value. We extracted the TWRs in the three-dimensional geometry of the aneurysm and the other region of the aneurysm wall was defined as non-TWRs. The mean values of each parameter in TWRs and non-TWRs were calculated and compared between the two groups using Mann-Whitney’s U test. As a result, the mean PD of all cases was 0.0688 in TWRs and -0.0278 in non-TWRs, the mean WSS was 0.0510 in TWRs and 0.0439 in non-TWRs, and the mean WSSD was 0.0162 in TWRs and 0.0014 in non-TWRs, respectively. These three parameters were statistically significantly higher for TWRs than non-TWRs (P<0.05). The higher PD of the TWRs is thought to be caused by the vertical stress on the aneurysm wall due to flow impingement on the aneurysm wall. We considered that high friction against the aneurysm wall is related to high WSS in the TWRs, resulting in a decrease the number of endothelial cells in the wall and the thin-wall. The WSSD was high in the TWRs because the tensile forces on the aneurysm wall may lead to the thin-wall. In conclusion, our results suggested that these hemodynamic parameters are related to TWRs. By conducting a CFD analysis on each patient and examining the areas with high PD, WSS, and WSSD values, it may be possible to identify TWRs before treatment.

Introduction

Image-based computational models of cardiac electromechanics (EM) are a powerful tool to understand the
mechanisms underlying physiological and pathological conditions in cardiac function and to improve diagnosis and therapy planning. To realize such advanced applications methodological key challenges must be addressed. First, enhanced computational efficiency and robustness is crucial to facilitate model personalization and the simulation of prolonged observation periods under a broad range of conditions. Second, physiological completeness encompassing therapy-relevant mechanisms is needed to endow models with predictive capabilities beyond the mere replication of observations.​​

Methods

In this talk, we report on a universal cardiac EM modeling framework that builds on a flexible method for
coupling a 3D model of cardiac EM to the physiologically comprehensive 0D CircAdapt model representing
closed-loop circulation. Additionally, we present recent advances in EM cardiac model personalization. In
particular, we focus on the identification of passive cardiac properties. Here, we present a novel methodology to simultaneously perform an automated identification of in-vivo passive mechanical properties and an estimation of the unloaded reference configuration.

Results

We report on the efficiency, robustness, and accuracy of the numerical scheme and solver implementation and show the model’s ability to replicate physiological behaviors by simulating the transient response to alterations in loading conditions and contractility, as induced by experimental protocols used for assessing systolic and diastolic ventricular properties. Further, we demonstrate the applicability of the framework to several clinically relevant problems.

Conclusion

The mechanistic completeness and efficiency of the framework renders advanced EM modeling applications feasible. The model facilitates the efficient and robust exploration of parameter spaces over prolonged observation periods which is pivotal for personalizing models to closely match observations. Moreover, the model can be trusted to provide predictions of the acute transient response to interventions or therapies altering loading conditions and contractility.

Coronary artery disease is a chronic disease in which the blood vessels are blocked, as a result of the growth of atherosclerotic plaques inside the arteries. To treat atherosclerosis, different clinical treatments are employed, including percutaneous coronary intervention (PCI). PCI includes the insertion of a stent, a metallic scaffold, its placement and expansion inside the diseased artery and the restoration of the blood flow. For a stent to be released to the market, its performance in terms of safety and efficacy is initially tested and evaluated, in in vitro, animal studies and clinical trials. However, this process is time consuming and costly, and there are certain limitations, such as the onset of severe complications and adverse events during clinical studies, which can be addressed by the exploitation of the very promising in-silico approaches. The fast-emerging technology of in-silico technologies is considered as an effective means for providing answers to “What if” questions. Among the available computational methods, the application of finite element method has proven to be the “method of choice” in replicating the biomechanical response of stents, by modelling their deployment in diseased, idealized or realistic, arteries. Structural simulations of stent implantation allow the evaluation of quantities of interest which could be difficult or even impossible to be observed through other means. By using models of coronary arteries and different stents or stenting techniques, the analysis of several geometrical (e.g., size of atherosclerotic plaques, malapposition, stent cell size, minimum lumen area, etc.) and mechanical quantities (type of atherosclerotic plaques, stent stress/strain, arterial stress/strain, etc.) can be performed. Such models provide useful insights into several aspects of stent design, towards of a final stent configuration for optimal post-implantation outcomes. In this work, a finite element analysis is performed to simulate the deployment of a commercially available stent inside a reconstructed patient specific artery. In more detail, for the reconstruction of the arterial wall and the underlying atherosclerotic plaque components, Optical Coherence Tomography and Invasive Coronary Angiography imaging are used, and an in-house 3D reconstruction software tool is employed. The next step involves the design of the 3D stent and its positioning in the stenosed artery. Then by the application of appropriate boundary conditions in the arterial-stent model and the implementation of a pressure driven approach, the in-silico expansion of the stent is achieved. For the finite model, representative material properties and material models are used. This study focuses on investigating the effect of the material properties of the underlying arterial wall and atherosclerotic plaques (fibrous, lipid, calcified) for providing guidance in appropriate material selection in stent modeling approaches towards improving stent design and performance. In more detail, we investigate the stresses, strains and deformations caused in the scaffold and the arterial wall considering the calcified plaques stiffness variance.

Among intermediate filaments populating cells, lamin filaments are of particular interest due to their universal role in the cell nucleus. Structural disruptions and mechanical alterations occurring in lamin mutations, commonly named laminopathies, result in increased cellular apoptosis and necrosis, leading to severe clinical symptoms. However, it remains quite unclear what determines the tensile stiffness and strength of alpha-helical rod domains in lamins, which are the elementary building blocks of the nuclear envelope. This study is concerned with the computation of these effective mechanical properties. Molecular dynamics simulations are used at the microscopic scale to predict the atomistic force fields of lamin fragments subjected to different axial stretches. The axial force of an equivalent continuous rod is deduced with the principle of virtual work. Using this methodology, we eventually predict the stiffness and strength of the rod for different configurations that are relevant for lamins. Results show that the tensile stiffness and strength of lamins are determined by non-local electrostatic interactions. Nonlocal interactions refer to points along the rod that are physically close in space despite not being parametrically close. Non-local electrostatic interactions result from heterogeneous distributions of electric charges, which are responsible for pretensions in lamin rods and play a major role in their tensile response. Our computational models show that mutations changing the charge of a single residue can significantly alter the individual mechanical behaviour of the lamin protein and hamper the ability of the nuclear envelop to withstand mechanical stresses. Such results pave the way towards the use of computational medicine to advance the treatments of laminopathies.

Spider silk is an important model system of a fully sustainable, natural biopolymer featuring multi-functionality and outstanding mechanical performance. Typical for a biomaterial, these properties are owed to a hierarchical structure spanning many decades from molecular to macroscopic scales. Because of its attractive characteristics, this material has been intensely studied in the past decades; nonetheless revealing its structure has proven surprisingly difficult. We combined atomic force microscopy (AFM), vibrational spectroscopy and other techniques to develop the most detailed structural model for a spider silk to date, using the recluse spider as a model system. Mechanical and nano-mechanical techniques were used to reveal macroscopic properties of the spider’s web and relate them to interactions between the constituents of the material, down to the level of specific bonds within the constituting amino acids.

At the centimeter scale, the spider web is “looped” into a 1-D meta-structure to increase toughness through a unique “strain cycling” mechanism. This meta-structure is held in place through a self-strengthening adhesive mechanism that avoids the weakness usually inherent to peeling/delamination failure. This inspired the design of a novel class of meta-materials with tunable mechanical properties. We further discovered that the ribbon silk of the recluse spider consists entirely of uniform, parallel nanofibrils about 10 nm in diameter and many microns long. The nanofibrils are relatively weakly adhering to each other, which gives the silk ribbon highly anisotropic mechanical properties. We developed a method combining AFM indentation with finite element analysis (FEA) to provide extensive characterization of the ribbon’s nanomechanical properties, including anisotropy of local stiffness and strength, as well as plastic deformation. The smallest length scale — at the level of spidroin molecules — was probed using vibrational spectroscopies: six different protein secondary structures were identified in silk nanofibrils for the first time, and their exact composition and angular orientation was determined. The level of detail revealed in our experimental studies will allow for computational models and analyses of this fascinating material in unprecedented detail. Multi-scale models can directly be validated and calibrated using experimental data at different length scales, with the prospect of developing a much better understanding of the material than before.

We further devised a disassembly process for silk fibers, which allows us to “reverse-engineer” silk fibers into their nanofibrillar components. This top-down approach allowed us to mass-produce native spider silk nanofibrils for the first time, with many applications and opportunities to further study nanofibrils. In addition, we developed a bottom-up approach where we visualize the molecular-scale assembly process of silk molecules into nanofibrils. In addition to structure, these methods also reveal details about the formation process of nanofibrils. Thus, we believe that these strategies will pave the way for synthesis of high-performance fibers and materials inspired by spider silk.

Breast cancer and prostate cancer tend to spread to the bone. Once they reach the bone site, cancer cells colonize and start to grow, which can have a negative impact on the bone. We have developed the first in vitro testbeds for studying breast cancer and prostate cancer metastasis to bone. Our research has shown that cancer cells undergo significant morphological changes as they cluster together and form tumors at the bone metastases site. These changes are related to the cancer phenotype and can dramatically alter the mechanical properties of the cancer cells extracted from the tumors. Eukaryotic cells consist of the cytoskeleton responsible for the cell's shape, internal organization, and mechanical rigidity. The cytoskeleton consists of actin microfilaments, microtubules, and intermediate filaments. Our studies on nanomechanical testing of bone metastatic breast and prostate cancer tumors using in vitro models indicate a progressive softening of the cells with metastasis progression. Confocal imaging of cancer tumors during progression has revealed quantitatively and spatially significant actin dynamics. Experiments showed the reduction and reorganization of actin filaments during metastasis progression. This work describes steered molecular dynamics (SMD) simulations to evaluate actin dynamics. Molecular nuances such as conformational locks and nonbonded interactions at the inter-strand and intra-strand interfaces regulate F-actin dynamics. We have observed significant changes in gene expressions related to actin and actin depolymerization factor cofilin (ADF/cofilin) associated with these changes Actin-related genes are downregulated, and cofilin genes are upregulated during metastasis. Steered molecular dynamics simulations of actin and actin with ADF/cofilin describe the mechanisms of the resilience of actin molecules and depolymerization of actin by ADF/cofilin.

Additionally, tumorigenesis occurs through integrin activation, a critical step in initiating the colonization of cancer cells at the bone site. Integrin protein acts as a mechanotransducer that establishes mechanical reciprocity between the extracellular matrix (ECM) and cells at integrin-based adhesion sites. This protein plays a critical role in cell-ECM adhesion and cellular signaling. As new biomaterials are being developed for tissue engineering applications, understanding cellular adhesion to engineered surfaces mediated by integrins is critical. We conducted steered molecular dynamics (SMD) simulations to investigate the mechanical responses of integrin αvβ3 with and without ligand binding for tensile, bending, and torsional loading conditions. The ligand-binding integrin confirmed the integrin activation during equilibration by opening the hinge between βA and the hybrid domain. This activation of liganded αvβ3 integrin influenced the molecule's stiffness observed during tensile loading. Furthermore, we observed that the interface interaction between β-tail, hybrid, and epidermal growth factor domains altered integrin dynamics. The deformation of extended integrin models in the bending and unbending directions of integrin reveals the integrin molecule's stored folding energy and directionally dependent stiffnesses. Along with available experimental data, SMD simulation results were used to predict the mechanical properties of integrin and reveal underlying mechanisms of integrin-based adhesion on polymer clay nanocomposite-based biomaterials. The evolution of actin dynamics and integrin activation is critical to metastasis and constructing a mechanobiological cell model of metastasis.

Globally, an estimate of 3.5 billion people is affected by oral diseases. Oral diseases disproportionally affect the most vulnerable and disadvantage populations and remains among the most common noncommunicable diseases. The range of diseases include dental caries, periodontal disease, tooth loss, oral cancer among many others. Untreated dental caries is the most common global health condition, affecting permanent teeth of 2 billion people and primary tooth of more than 500 million children. An estimated 800 million resin composite, 100 million amalgam, and millions of glass ionomer cement restorations are placed each year and are one of the most prevalent medical interventions in the human body, not to mention the over five million implants placed in the United States each year. The cost for these therapies is immense. The combined complexity and prevalence of dental diseases requires well engineered and effective interventions to optimize patient outcomes. Our particular focus has been on exploring biomimetic approaches that can impart the key biological activity toward prevention of dental and oral diseases as well as restoration of oral health.

Biomolecules are essential in physiological processes, and smaller molecules such as peptides offer to mimic their function to be uniquely versatile molecular tools in coupling biologically instructive roles in material systems. We have been exploring peptide-based approaches including repairing mineralized tissues of teeth damaged due to caries, trauma or periodontal diseases and derivatizing dental materials with wide range of bioactivities including antimicrobial activity toward preventing infection. We identified several intrinsically disordered peptides and utilized the resulting sequence information to recapitulate the functional domains of native proteins. Several of these peptides demonstrated a rich conformational propensity that can be further tunable at the bio-hybrid material interfacial design. To address the different unmet health care need, we designed peptides with different properties ranging from self-assembly to mediating mineralization or catalytic function. By adapting experimental and computational approaches combined with transparent machine learning (ML) methods, we continue our search for the peptides that present desired function. In our approach, we expanded our predictions from single function to multifunctional peptides to control their properties at the complex materials to dental tissue interfaces. In this presentation, we will discuss our approach for next generation treatments that are composed of peptides, peptide-mimetics and peptide-polymer hybrids. We designed antimicrobial peptides (AMPs) and developed localized delivery approaches for their easy deployments on the sites. We developed machine learning approaches to classify antimicrobial peptides and enrich sequence domains for effective antimicrobial search with enhanced properties. By incorporating multi-pronged strategy to repair and protect exposed dentin and collagen, we combined different peptide, peptide-mimetic and peptide-polymer hybrid to treat different lesion types. Our peptide design and engineering approach targeting biohybrid interfaces provides an alternative delivery strategy to deploy peptides on the sites and increases their availability and preserving their efficacy. Peptide hybrids could offer multi-pronged and microinvasive strategies to overcome the limitations of current approaches as well as address the high prevalence of root caries and dental erosion in the aging population.

The expanded use of dental implants has led to an increase of peri-implant disease that shortens implant life and leads to failure. Peri-implant disease results from microbiota dysbiogenesis triggering in the host an immune inflammatory response that destroys soft and hard-tissue. At an incidence of 14.5%, and over 3 million implants placed and growing by 500,000/year, a reduced service life ending in failure will adversely impact public health and increase health care costs. Molecular studies on proteins expressed during tooth formation led to insights that can be translated to clinical care to improve health, such as development of an antimicrobial bifunctional peptide film to slow disease progression. Based on a high-affinity titanium binding peptide that anchors the anti-microbial peptide to the implant surface, greater than 98% coverage is achieved in <2 minutes even in the presence of contaminating protein, to produce a film durable to mechanical brushing that kills >90% of bacterial. Also developed was a short amelogenin peptide M59, that activates osteogenesis through the Wnt pathway. M59 can also be applied locally using the titanium binding peptide to deliver the M59 osteogenic signal where needed and on command. This non-surgical approach can improve oral health by delivering a simple to apply/reapply multiple times, using a water-delivered bifunctional peptide film serving to control microbial dysbiogenesis, reduce disease progression and induce bone formation.

According to the World Health Organization, 375,304 deaths and about 1.4M cases of prostate cancer reported in 2020. Majority of these deaths occur due to bone metastasis, or the transport of cancer from prostate to bone and subsequent skeletal failures. The highly diminished availability of bone metastasis prostate cancer samples and failure of animal models due to death preceding metastasis in animals necessitates the development of realistic prostate cancer in vitro models. A novel nanoclay-based tissue-engineered polymeric scaffold sequentially seeded with human mesenchymal cells and cancer cells recapitulates prostate cancer bone metastasis. We hypothesized that interstitial fluid flow acts as a driving force for migration of cancer cells in the vicinity of capillary pores. The fluid-flow provides unique biochemical cues for metastasis. Earlier we have reported the role of interstitial fluid in prostate cancer bone metastasis using a 3D in vitro bone metastasis testbed, with a perfusion bioreactor. However, understanding the metastatic cascade of cancer cells is essential to pave the way to discover therapies for metastatic cancers. In particular, the extravasation stage, when cancer cells start to invade into the secondary site resulting is metastatic tumors is critical. Extravasation comprises of first the transmigration across the capillary’s endothelium. Hence, we designed a bioreactor that enables interstitial fluid-flow for prostate cancer bone metastasis. We developed a novel 3D in vitro dynamic horizontal bioreactor integrated with transwell inserts, that recapitulates the in vivo microenvironment of cancer cells representing migration of cancer cells under interstitial fluid flow. Computational fluid dynamics studies, detailed biological characterizations including evaluation of gene expressions and imaging were conducted to build relationships between the fluid induced shear stresses and the progression of cancer metastasis. The computational fluid dynamics results indicate that 0.05 ml/min flow rate recapitulates the physiological condition. We evaluate the migration of the highly metastatic PC3 prostate cancer cells through the transwell insert under both dynamic and static culture conditions. While experimentally subjecting the cancer cells to the fluid derived stresses via the bioreactor experiments, we investigate the molecular mechanisms responsible for the migration of cancer cells under different culture conditions. This study demonstrates that adhesion proteins avb3 integrins play important roles in response to mechanical cues and act as mechanosensory agents that transport mechanical signals via the avb3-MMP 9 signaling axis to promote flow-induced motility of prostate cancer cells. Also, the interstitial fluid-flow does not alter the CXCR4 level, an important regulator of metastasis and invasiveness, but bone proximity upregulates CXCR4 levels enabling increased MMP-9 levels. Further, avb3 integrins and MMP-9 levels are upregulated by fluid-flow causing increased migration under fluid-flow. Overall, these studies describe the critical role of interstitial fluid-flow in prostate cancer metastasis.

Purpose

The meniscus is a critical component of a healthy knee joint and plays a pivotal role in preserving the knee homeostasis and biomechanics. Consequently, meniscal damage affects the knee equilibrium, progressively contributing to cartilage disruption up to the development of osteoarthritis (OA). Considering the meniscal poor-self healing potential, its repair still represents a challenge to orthopedic surgeons. To counteract the increasing demand for meniscus implants and allografts, innovative, and effective repair strategies are required and among vanguard approaches, 3D printing technologies seem to be intriguing and promising.

Material and Methods

To pave the way towards the development of 3D printed patient-specific menisci, we conducted an experimental study to, (i) create a virtual 3D reconstruction of the meniscus microstructure based on μCT scans; and (ii) build a STL model of the meniscus for 3D printing based on the scanned models.

The menisci used in this study were removed from TKA human knees. All tissue samples were harvested immediately after surgery, and a staining and freeze-drying protocol with Lugol’s solution was applied. μCT scans of the menisci were performed using a SCANCO μCT 50 (SCANCO Medical AG, Brütistellen, Switzerland) specimen scanner. A 3D voxel model was reconstructed from the data and converted to STL format using MIMICS (Materialise, Leuven, Belgium).

Results

A Lugol stained and freeze-dryed meniscus can be scanned via μCT. The μCT images clearly displayed microfibers in the meniscus with an isotropic resolution of 20.7 μm. The surface layer, lamellar layer, circumferential fibers, and radial fibers could be identified. Furthermore, a 3D STL-model of the meniscus was digitally built based on the μCT images, and the microscale fibers could be discriminated in the virtually reconstructed meniscus.

Conclusion

The microstructure of the meniscus can be visualized via μCT and can be reconstructed into STL format. The 3D STL model will be the data basis for generating a printable model of patient-specific menisci via 3D bioprinting.

Introduction

The knee joint is one of the most complex human joints and injuries like ligament rupture, cartilage or meniscus defects are common. Instabilities caused by these defects or suboptimal treatment lead to secondary complications. Patient-individual models that take into account the patient-specific knee joint motion are supposed to largely improve surgery planning and outcome. State-of-the-art finite element (FE) knee joint models are highly complex, computationally demanding and time-intensive.

Methods

We propose a simplified knee joint model based on position-based dynamics (PBD) for real-time prediction of patient-individual knee joint motion [1]. The model’s outcome is compared to the predicted patient-individual knee joint motion of a state-of-the-art FE model [2]. Both models are automatically generated from magnetic resonance (MR) images of a patient’s knee. Comparison is done using MR images of eleven healthy volunteers. The study was approved by the ethics committee of the Albert-Ludwigs University Freiburg (Nr. 91/19 – 210696, 19 August 2021) and all volunteers gave written informed consent prior to participation. We compared two motion sequences going from 0° knee flexion to 20° flexion and adding 1) internal or 2) external rotation torque. Both methods have previously been validated against reference MR measurements.

Results and Discussion

Important quality measures are the accuracy of the models as well as their computational performance, such that the clinical workflow is not interrupted and fast, physics-informed surgical planning can be performed.
A direct comparison of the FE and PBD model results shows similar motion predicted by both methods. Calculating one of the described motion sequences for one patient takes about 2 h using the FE approach. Applying the PBD model, one motion sequence can be calculated in 48 s on average.
Cartilage contact areas in the knee joint during motion are estimated and compared using both approaches. However, questions concerning the stresses and strains in ligaments, cartilages or the menisci can only be answered using the more complex FE model.

Conclusions

From the above observations, we see that the high modelling accuracy of FE knee joint models is dispensable in our use case of predicting knee joint motion. The PBD model predicts similar motion while being approximately 150-times faster. Thus, the PBD model is even capable of real-time modifications, e.g., repositioning of a planned anterior cruciate ligament transplant. Summarizing, the PBD model yields a good trade-off between accuracy and time constraints.

Acknowledgements

This work has been supported by the German Ministry of Education and Research (BMBF) in the call “Individualisierte Medizintechnik” under contract number 13GW0277. We would like to thank Ingmar Ludwig, Elin Wefer, Elham Taghizadeh, Sebastian Bendak, Jonas Buchholtz, and Hagen Schmal for their support of this work.

References

[1] I. Ludwig, E. Taghizadeh, K. Izadpanah, T. Lange, J. Georgii, “Patient-specific Modelling and Simulation of Knee Joint Motion using Position-Based Dynamics”, in Proc. of CARS, 2020.

[2] E. Theilen et al., "Validation of a finite element simulation for predicting individual knee joint kinematics," in IEEE Open Journal of Engineering in Medicine and Biology, doi: 10.1109/OJEMB.2023.3258362.

Children with limb deformities often seek paediatric orthopaedic consultation because of uneven growth of the physis or growth plate of the bone structure. This process can be influenced by mechanical stimulation. Surgical implantation of devices such as staples can restrict growth in specific areas of the affected cartilage to correct the asymmetry. Although these techniques temporarily block the physis, they are less effective in certain deformities and can cause complications. The development of computational growth models will improve these implant-guided growth techniques. Therefore, our aim is to develop a computational tool to predict epiphyseal growth and support clinical strategies to correct abnormal growth of immature bone.

The model developed is based on the mathematical model of Garzón-Alvarado et al.¹, which is a biologically derived model that considers both the mechanical and cellular activity of the growth plate. The growth rate is described by a tensor, which is the sum of the contribution of both proliferating and hypertrophying chondrocytes to the increase in height in the growth direction. Each term has a biological growth component and another component describing the effect of increased hydrostatic and deviatoric stresses associated with growth in the respective zones of the physis. This model has a significant advantage in that it describes the pattern of chondrocytes in each zone using a tensor. We have taken the parameter values from the literature and the stress-related values are approximations. In addition, we have introduced a parameter in the hypertrophic zone that is proportional to the relative growth of the proliferative zone. This parameter allows to block the physis under high stresses, as occurs when a staple is inserted.

We compared the growth generated by our model with the results published by Narváez et al.² for the case of proximal tibia growth in rats over a 23-day period. Similar growth values were obtained in the simulation of free growth, i.e.where there is no increase in stress due to external loads. Our model also showed a marked tendency for the growth rate to decrease with the application of compressive loads as well. We also simulated Narvaez's experiment in which the growth plate was subjected to traction during the 23-day period. Our results differed from the original experiment in that they produced less growth and varied the geometry of the plate. This could be explained by the hypothesis that the growth is mainly due to the ability of the cells to proliferate and then hypertrophy. Therefore, if distraction does not have a positive effect on proliferative activity as some experiments have shown and may even decrease it^3,4, we are inclined to think that our results are physiologically relevant.

References

1. Garzón-Alvarado et al., J. Mech. Med. Biol 11(5): 1213-1240 (2011)
2. Narváez-Tovar et al., Theor Biol Med Model 9, 41 (2012)
3. Alberty et al., Acta Orthop Scand 64(4):449-55 (1993)
4. Apte et al., J Bone Joint Surg Br 76(5):837-43 (1994)

Acknowledgements

Spanish Government and Health Institute Carlos III (PI20/00293); Children Hospital Sant Joan de Déu; Department of Communication and Information Technologies, Universitat Pompeu Fabra.

The human shoulder joint combines mobility and stability in a unique musculoskeletal system. The anatomical structure of the glenohumeral joint allows for an extensive range of motion, while passive and active soft tissues ensure the joint’s integrity through static and dynamic mechanisms. In this context, muscles, especially the rotator cuff and the deltoid, play an essential role. On the one hand, muscles actively stabilize the otherwise extremely unstable glenohumeral joint through the mechanisms of concavity compression and scapulohumeral balance. On the other hand, muscles act as torque generators and enable complex movement patterns through their sophisticated interplay. Maintaining this delicate balance between mobility and stability is essential for proper shoulder function, yet it is easily disrupted by injury or pathological conditions. Despite the high incidence of shoulder disorders in clinical practice, understanding of the underlying biomechanics remains limited, posing major challenges for medical assessment and treatment.
Computational musculoskeletal models offer great potential for biomechanical studies of the shoulder’s physiology, investigations of pathological conditions, objective predictions and evaluations of (patient-specific) treatments, and the development of rehabilitation equipment for physical therapy. While numerous reduced-dimensional multi-body models exist, research on comprehensive continuum-based finite element models remains limited. However, three-dimensional interactions between the joint components, such as contact and sliding mechanisms, are central to the shoulder’s physiology. In contrast to multi-body models, continuum-mechanical models can represent such volumetric effects, account for complex muscle fiber and tendon arrangements, and model sophisticated constitutive behavior. For biomechanical studies of the shoulder, they are thus particularly relevant.
Considering their role as active joint stabilizers and force generators, skeletal muscles deserve special attention regarding their material description. Passive skeletal muscle is - according to its histological composition - commonly modeled as a transverse isotropic composite of unidirectionally oriented fibers connected by extracellular tissue. Active contractile effects are incorporated through active-stress or (generalized) active-strain approaches. Since selecting an appropriate constitutive model is crucial for reliable predictions, the question arises of which material is best suitable for characterizing the shoulder’s skeletal muscles.
In this contribution, we contrast three hyperelastic formulations considering mathematical, computational, and physiological aspects: an active-stress, an active-strain, and a generalized active-strain approach. To establish a basis for comparison, we fit the material parameters to a common set of experimentally obtained stress-strain data. As one load case is generally insufficient to determine the material response uniquely, we consider multiple active and passive loading conditions. We discuss the concepts of modeling active material behavior from a mathematical and physiological perspective, address analytical and numerical problems arising from the mathematical formulations, and analyze the included biophysical principles of force generation in terms of physiological correctness and relevance considering the modeling of the human shoulder. Conclusively, we present a constitutive model combining the studied materials’ most promising and relevant properties. By the example of a fusiform muscle geometry, we investigate force generation, deformation, and kinematics during active isometric and free contractions. Eventually, we demonstrate the applicability of the material formulations in simulations of a comprehensive continuum-mechanical model of the human shoulder.

Surgical resection of chest wall tumours or large defects due to thorax trauma with rib fractures may lead to a loss of ribcage stability and require reconstruction to allow for physical thorax functioning. When titanium implants are used especially for larger, lateral defects, they tend to fail due to implant fracture, peri-implant rib fracture or screw loosening. Implant failures are mainly caused by the specific mechanical requirements for chest-wall reconstruction which must mimic the physiological properties of the ribcage and which are not yet met by available implants.

Rib implants must show some important characteristics: On the one hand, there are a variety of quasi-static, dynamic and fatigue loads that they must be able to withstand permanently. Quasi-static loads can be induced by body posture, e.g., lying on the side, whereas dynamic loads are mostly induced by impacts. The most important fatigue load is normal breathing with approximately 8.4 million cycles per year. On the other hand, these implants must have an appropriate stiffness and strength to enable all daily movements and at the same time protect the vital inner organs. To meet these requirements, it is essential to understand the biomechanics of the thorax. For this purpose, a full thorax FEM model was developed, comprising all biomechanical relevant structures, e.g., ribs, sternum, costal cartilage, vertebral column, costo-vertebral joints, vertebral disks, passive intercostal muscles, subcutis and skin.

The simulation model was assembled in a stepwise approach. First a chest CT scan of a fresh, unembalmed cadaver in the supine position was made to reconstruct the anatomical structures of bones and cartilage. Additional CTs in different positions as well as stiffness measurements on several anatomical structures and levels were used to define mechanical properties of the FE model and for calibration and verification. Various activities such as ventilation, breathing, resuscitation, lying on the side or coughing were simulated on this verified FE model and the deformations and stresses were evaluated. Several simulations were carried out with different defects on the rib levels 5 to 9 and with corresponding implants according to the current state of the art. the critical activities that lead to damage of the implants could be identified using the prescribed procedure. Based on these findings, an algorithm for determining the implant dimensions for different alternative metallic and non-metallic materials (PEEK, …) was developed to achieve the required stiffness and strength of 3D-printed rib implants. This paves the way to optimized, patient-specific and intraoperatively 3D-printed rib implants.

Background: Cephalomedullary nailing is frequently used to treat per- and subtrochanteric fractures of the proximal femur. After fracture healing, patients sometimes request nail removal due to persistent pain or irritation. However, removing the nail leaves a large void in the bone, which poses a considerable risk of re-fracture at the femoral neck. Pre-operative prediction of fracture risk would help to make an informed decision about nail removal and to estimate the required post-operative care. This study investigated whether patient-specific finite element (FE) models created from pre-operative CT scans can predict femoral bone strength after nail removal. Experimental data of femora after nail removal were used to evaluate the accuracy of the models.

Methods: Ten femoral bones of human body donors who were treated with a cephalomedullary nail during their lifetime due to a per- or subtrochanteric fracture were obtained from the Medical Bio-/Implantbank Vienna. CT scans (0.4x0.4x0.6 mm³ voxel size) were taken prior to nail removal using a dual energy protocol and an iterative metal artefact reduction algorithm. The bones were cut to 50 % length, embedded, and mounted to a material testing machine to simulate loading in stance. The load was increased monotonically until failure and the maximum force was recorded. The experiments were replicated using patient-specific nonlinear voxel-based FE models. The models were created by virtually removing the implant from the pre-operative CT image, aligning the bones in agreement with the experiment, coarsening the image to 3 mm voxel size and converting each voxel to a linear hexahedral element. Due to remaining metal artefacts from the distal locking screws, the FE models were cut to the proximal region above the distal locking screw. A density-dependent, isotropic, elastic-damage material was assigned to each element and the models were loaded until failure in analogy to the experiments. The maximum force predicted by the models was then compared to the experimentally measured maximum force using linear regression analysis.

Results: Experimental femoral bone strength after nail removal ranged from 611 to 2851 N and FE-predicted strength ranged from 390 to 1873 N. FE model predictions and experimental measurements were well correlated (R²=0.78, p<0.001), but the models underestimated the experimental measurements (experimental mean: 1837±598 N, FE mean: 1127±425 N).

Conclusions: The FE models were able to predict the strength of femoral bones after cephalomedullary nail removal pre-operatively with good correlation to experimental measurements. This shows that voxel-based FE models can predict bone strength despite the presence of a metal implant in the CT scan and the highly irregular structure of the previously fractured and healed bones. Thus, FE models may be a useful tool to support clinical decisions on nail removal in the future.

Introduction: Nonlinear micro-finite element (μFE) models represent a powerful tool to predict elastic and yield properties as well as damage onset of bone across length scales ¹. Apparent mechanical properties and damage predictions have been validated against experimental measurements at the macroscale ^2,3. Recently, validation of local properties and damage onset have been achieved using digital volume correlation (DVC) based on micro-computed tomography (μCT) images in trabecular bone ⁴. However, validation of damage predictions in cortical bone remains missing due to the limited μCT studies capturing the deformation of cortical bone at high spatial resolution. Here, we aim to (i) use experimentally measured 3D displacement fields to validate nonlinear μFE models of cortical bone; (ii) use those models to investigate damage emergence and propagation based on synchrotron-radiation (SR)-μCT imaging.

Materials and methods: We performed time-resolved in situ SR-μCT compression testing in bovine cortical bone specimens in beamline I13-2 at Diamond Light Source. Specimens were loaded while SR-μCT images (6.5 μm voxel size) simultaneously acquired during compression up to apparent failure. We used DVC to obtain full-field displacement fields in the specimens ⁵. We generated μFE models using the unloaded SR-μCT image coarsened to 26 μm voxel size, and applied DVC displacement fields as boundary conditions ⁴. An elastic-viscoplastic damage model ¹ featuring an isotropic Drucker-Prager yield surface ⁶ was used (UMAT, Abaqus v6 R2018), with an isotropic Young’s modulus of 22.8 GPa ⁷ and a Poisson’s ratio of 0.3. Correlations between DVC measurements and μFE predictions in the apparent elastic regime and prior to failure were investigated using concordance correlation coefficient (r_c). We validated damage predictions against cracks location from the SR-μCT images post-failure.

Results: Displacements predictions using nonlinear μFE models outperformed linear μFE models both at the apparent elastic (r_c,linear ≥ 0.71, r_c,nonlinear ≥0.85) and plastic (r_c,linear ≥ 0.74, r_c,nonlinear ≥0.87) regions. Nonlinear μFE models damage predictions seem to initiate next to vascular porosity and correlated to regions that displayed significant cracks post-failure in the SR-μCT images.

Discussion: Our results demonstrate the ability of nonlinear μFE models to accurately predict displacements and to capture damage location in cortical bone tissue. They can, thus, be explored to investigate the mechanisms of bone failure in relation to structural and material changes due to aging or disease, enabling the development of treatment strategies that prevent bone fracture.

References: 1 Schwiedrzik, J. J. et al. Biomech. Model. Mechanobiol. 12, 201–213 (2013)e. 2 Hosseini, H. S. et al. J. Mech. Behav. Biomed. Mater. 15, 93–102 (2012). 3 Dall’Ara, E. et al. Bone 52, 27–38 (2013). 4 Peña Fernández, M. et al. J. Mech. Behav. Biomed. Mater. 132, 105303 (2022). 5 Peña Fernández, M. et al. Acta Biomater. 131, 424–439 (2021). 6 Schwiedrzik, J. J. et al. Biomech. Model. Mechanobiol. 12, 1155–1168 (2013). 7 Schwiedrzik, J. et al. Acta Biomater. 60, 302–314 (2017)

Acknowledgements: Leverhulme Trust RPG-2020-215

Collagens are the primary constituents at the lowest hierarchy of biological tissues (Orgel et al., 2006). Abnormality in the stiffness of collagen fibrils may result in severe cardiovascular diseases, such as atherosclerotic plaques (Akyildiz et al., 2017). In this context it is an open quest to describe the biomechanical properties of collagen fibrils by a constitutive model.

Hyperelastic models are a powerful approach to depict the nonlinear mechanical response of biological materials (Holzapfel & Ogden, 2020). However, hyperelasticity fails to describe non-affine transformation and inherent rate-dependency observed in soft tissues (Morin et al., 2021). In this study, a data-driven hypoelastic model for individual collagen fibrils is attempted to be built within the framework of continuum theory, which allows in-depth discussion on fiber-rearrangement and rate-dependent nonlinear elastic behavior. The experimental data are provided by a previous study of atomic force microscopy (AFM) uniaxial tensile tests on individual collagen fibrils from rat tail tendons (RTT) in low strain regimes (Andriotis et al., 2018). Owing to the large aspect ratio of the collagen fibrils and the rotational-free tensile tests, one can infer the assumptions on linear displacement field and zero spin rate for further stress- and strain-rate calculation. Hence, the resulting hypoelasticity tensor is defined by the Cauchy stress-rate and strain-rate field, which can be derived from the force-displacement data of AFM tensile tests. With the absent knowledge on the cross-sectional area of the deformed collagen fibrils, at the current stage this research focuses on the analysis of the measured displacement, strain and force rate as well as the computation on the one-dimensional extensional stiffness.

Results show that the measured displacement rates of the collagen fibrils increase with the level of displacement, and deviate from the displacement rates originally set in the experiment, being up to 1.72-fold higher. Also at low strain, nonlinear extensional stiffness of collagen fibrils is observed. As both the force rate and the extensional stiffness reach a plateau at the strain level of approximately 2.6 from each experiment, it is suggested that the tensile response of the collagen fibrils is reaching the maximum extensional stiffness of phase 1 deformation. The force-displacement profiles with respect to different set displacement rates may indicate viscoelastic behavior of collagen fibrils. Such that, the viscoelastic contribution can be supplemented into the developing hypoelastic model with the addition of a dissipation function (Rajagopal & Srinivasa, 2011). Subsequently, this research will proceed in the estimation of the three-dimensional hypoelasticity tensor with the inclusion of collagen fibril ultrastructure.

Growth in soft tissues is a continuous process of mass deposition/reabsorption to maintain a homeostatic mechanical state. While the classical g rowth models aim at quantifying the evolution of the tissue mechanical response, this contribution proposes an extension towards the realm of truly multiscale models so as to capture local strains and stresses during growth.

In more details, a representative volume element (RVE) of the soft tissue is considered which is made of growing spherical inclusions embedded in a soft matrix, both phases displaying a hypoelastic constitutive behavior [1]. The RVE is subjected to a uniform macroscopic strain rate over a given time period. Then, we translate the (Kröner-Lee) multiplicative decomposition of the deformation gradient to the microscopic deformation gradient averaged over the volume of each growing phase. According to [2], the microscopic strain in each phase is then computed as a linear function of both the imposed macroscopic load and the inelastic growth-related (eigen)strain. The latter evolves according to a properly-defined evolution law, between the growth-related deformation gradient and the mass production rate. In our case, we consider the growth to be isotropic, i.e. proportional to the identity tensor. Besides, the net mass production rate is set to be proportional to the difference between current and target homeostatic stresses [3]. Finally, we obtain the macroscopic consequences on the tissue (e.g. in terms of tissue growth and deformation) from average rules. Different case studies are investigated, namely the deposition of a new phase within the RVE and the progressive degradation of an existing phase.

This study shows that our model is able to reproduce the trends of growth-induced changes in terms of phase volume fraction and local deformation fields. Such an evaluation of local stress and strain fields is of particular interest to evaluate the mechanical environment sensed by populations of mechanosensitive cells, contained in the tissue. Therefore, accounting for mechanobiology will be a necessary future work to capture the mechanisms of growth in the context of aneurysms for instance and enable the calibration of model parameters.

[1] Morin, C., Hellmich, C., Nejim, Z., and Avril, S., ”Fiber Rearrangement and Matrix Compression in Soft Tissues: Multiscale Hypoelasticity and Application to Tendon”, Front. Bioeng. Biotechnol., 9:725047 (2021).

[2] Pichler B, Hellmich C. ”Estimation of influence tensors for eigenstressed multiphase elastic media with nonaligned inclusion phases of arbitrary ellipsoidal shape”. J Eng Mech. 136(8):1043–1053 (2010).

[3] Braeu, F. A., Seitz, A., Aydin, R. C., and Cyron, C. J., ”Homogenized constrained mixture models for anisotropic volumetric growth and remodeling”, Biomech. Model. Mechanobiol., 16(3), 889-906 (2017).

Early detection and mechanical characterization of atherosclerosis plaques can aid in determining a patient’s vulnerability to acute coronary syndromes. Currently, intravascular ultrasound (IVUS) is the most common medical imaging methodology for the detection of atherosclerotic plaques. In combination with new techniques, it provides an approximated field of strain or displacement of the coronary wall. In this work, we present two methodologies to characterize the mechanical properties of atherosclerotic tissues. Both methods use the approximated radial strains in an optimization process to estimate the material parameters of the tissues. The first methodology is based on a classical estimation of Young’s modulus of the tissues, while the second one focuses on estimating the hyperelastic material properties and the geometry of the unpressurized plaque. IVUS data were simulated by in silico models to implement and validate the methodologies. We analyze different materials and real geometries in the finite element (FE) models. We add a signal-to-noise ratio of 20dB to the strain field in order to simulate the intrinsic noise present in IVUS data. The segmentation is based on the representation of strain variables and a Watershed process to extract the different tissues [1]. Then, an inverse FE analysis is performed to estimate the mechanical behavior of the tissues. Both methodologies use a pattern search algorithm for the optimization process and take a maximum of four hours. The first methodology estimates the relative stiffness of the tissues at certain blood pressure [2], whereas the other uses a Pull-Back algorithm [3] to recover an estimated unpressurized geometry to obtain the hyperplastic material response. The second method significantly improves over the first, resulting in a 13% error reduction in the radial strain map of the optimization process. In addition, the second method provides a more detailed explanation of the mechanical response throughout the entire physiological cardiac cycle. Although the estimation of fibrotic tissue materials is accurate, the properties of calcifications and lipids are less accurate. These results indicate that the fibrotic tissues, rather than the lipid core or calcification, mainly influence the strains observed in atherosclerotic plaques. In the second methodology, we develop a new technique that successfully estimates the hyperelastic material parameters for atherosclerotic tissues and provides an unpressurized geometry. The findings enable the estimation of the stress field over the plaque, providing valuable insights into its vulnerability to rupture. Although this work is purely theoretical, a new pipeline is defined and its preliminary results are promising.

Acknowledgments

This work was supported by the Spanish Ministry of Science and Technology through research project PID2019-107517RB-I00, the regional Government of Aragón support for the funding of the research project T24-20R, and grant CUS/581/2020. CIBER Actions are financed by the Instituto de Salud Carlos III with assistance from the European Regional Development Fund.

References

[1] Latorre ÁT. et al. Mathematics. 10 (21), 4020, 2022.

[2] Le Floc’h S. et al. IEEE Trans Med Imaging. 28 (7), 1126–1137, 2009.

[3] Raghavan ML. et al. Ann Biomed Eng. 34 (9), 1414–1419, 2006.

Atherosclerosis is the most common cardiovascular disease, but also the most common cause of other cardiovascular diseases. Carotid arteries are the primary source of oxygenated blood to the brain, and their narrowing (stenosis) by atherosclerotic plaque can limit blood flow in the brain. Plaque rupture can lead to a complete blockage of blood flow, which in turn can result in a life-threatening stroke or myocardial infarction. There are two options for surgical treatment: open surgery or endovascular treatment (angioplasty), which is often combined with a stent implementation. According to the method of deployment, stents are divided into balloon-expanding (e.g., metal stents) and self-expanding ones (e.g., made of shape memory alloys). During the installation of balloon-expandable stents, unnecessarily high loads are applied to the arterial wall, which can cause injuries.

With the development of new materials, new stents are continuously being developed with the aim of reducing complications caused by their deployment, such as in-stent thrombosis or restenosis. Polymer stents are known to have lower resistance to radial compression compared to traditionally used metal stents, however, they are increasingly used in the production of medical implants because they are associated with easier recovery, faster and cheaper production, biodegradability of materials, and reduction the invasiveness of the deployment procedure.

As shown by various researchers, but also this work, the geometry of the stent can directly influence the reduction of the possibility of complications, and thus increase its effectiveness. However, optimizing the stent construction is an extremely complex process.

The primary purpose of this work is the efficiency analysis of bionic-inspired stent models made of shape memory polymer in comparison to a conventional stent with reference to residual stenosis. Since the real artery wall behaviour is anisotropic endowed by the arrangement of collagen fibres, the artery wall was modelled as an anisotropic hyperelastic material described by the Holzapfel-Gasser-Ogden (HGO) model and isotropic hyperelastic material described by the Neo-Hooke model was used for atherosclerotic plaque and stent. User material subroutine (UMAT) was used for the implementation of HGO model with elastin and collagen pre-stretch in order to define maximally realistic loading conditions of the artery.

Efficiency analysis was achieved consecutively by defining a conventional stent model, proposal of bionic-inspired stent models and finally by comparing the numerically obtained results of the observed parameters in the case of stiff hypocellular atherosclerotic plaque whose composition and mechanical properties generate the most threatening case of the disease, smallest degree of stenosis reduction and greatest artery wall stress. It was determined that the model inspired by the honeycomb structure provides more favourable conditions for stent implantation in comparison to the conventional and other proposed bionic-inspired models.

Acknowledgement:

This work was supported, in part, by a grant from the Croatian Science Foundation (IP-2020-02-4016, PI: Ž. Tuković).

Clinical decision making for abdominal aortic aneurysm (AAA) treatment is based on the measurement of the maximum diameter from medical images (Wanhainen et al., 2019). However, the diameter alone has shown a limited ability to correctly predict the risk of rupture, which in turn can negatively affect the patient’s survival rate (Vorp, 2007). Biomechanics-based indexes have proven better predictive ability, but they have not yet found practical application in clinical context (Polzer et al., 2020). In vivo non-invasive estimations of AAA tissue properties can help to assess the AAA vulnerability. Methods based on time-resolved ultrasound (US) were proposed, US being the most commonly employed imaging technique for AAA monitoring (Van Disseldorp et al., 2016). However, they should take into account the effect of surrounding tissues (Petterson et al., 2019). It was also noticed that US scanning procedure itself introduces an external force applied to the AAA via the transducer push, causing deformations in the aortic wall (Ghulam et al., 2022). In this work, we propose to correct the AAA material properties estimation by means of a multi-factorial correction coefficient, taking into account the geometry, the surrounding tissues and the probe pressure thanks to finite element modelling. Three patient specific geometries of AAA and surrounding tissues were obtained from clinical images via segmentation. Specifically, the soft tissues surrounding the spine up to the patient surface were simplified as a homogeneous material. The spine geometry was subtracted from the soft tissues and it was modelled as a boundary condition, fixating the internal surface created by the cut. Through this approach, both the effect of the hard and soft surrounding tissues were considered. The AAA wall was modelled as orthotropic linear elastic material as in Perrin et al., 2015, while the surrounding tissues were modelled as hyperelastic neo-hookean materials as in Petterson et al., 2019, and the value of shear stress was obtained by calibrating the simulation according to the data presented in Ghulam et al., 2022. The inner wall of the AAA was subjected to patient-specific blood pressure and the outer wall was tied to the surrounding tissues. We also applied a localized uniform firm pressure of 13.85 kPa (value from Ghulam et al., 2022) on the outer surface of the patient in the antero-posterior direction. The stiffness of the AAA wall was then calculated as the ratio between the hoop stress, calculated with the Laplace relationship, and the hoop strain. The effect of geometry was quantified by comparing the results for the pressurized patient specific AAA geometry to an idealized case, i.e. a pressurized thin walled cylinder, obtaining the first correction factor nu_G. Progressively adding complexity, we obtained two more correction factors nu_ST and nu_PP respectively correcting for the presence of surrounding tissues and probe pressure. We found the correction factors to be consistent across patients: nu_G = 0.70 (0.042), nu_ST = 2.01(0.002), nu_PP = 1.14 (0.025). We conclude that the presented method is a promising and simple approach to account for the different factors affecting US-based biomechanical estimations.

This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 859836.

The segmentation of computed tomography (CT) scans allows to define patient-specific computational models for the development of targeted tissue engineering applications. However, the segmentation of these images does not provide an accurate definition of the soft tissues. For cervical spine models, this problem is observed when rebuilding intervertebral discs (IVD), facet joints and ligaments, requiring adjustments on the geometry of these structures. In the literature, there is no consensual definition of the parameters that characterise these structures, as they provide different results when compared to the benchmark experimental data. Also, recent works use non-linear material properties to achieve a better physiological cervical spine mobility. The aim of this study was to define the role of the components of the cervical spine, on a linear finite element (FE) model, and how can a combination of modelling parameters be adapted to better represent cervical spine mobility and assist on enhancing new tissue engineering strategies. Adjustments on IVD geometry, ligaments insertion and material properties were taken to evaluate if a FE model could replicate experimental data mobility without requiring non-linear material properties.

Three subject-specific models were segmented and converted to FE models in ScanIP (Synopsys, Inc.,USA). The IVD design was based on vertebrae boundaries and experimental data volumes. The segmentation of the IVD was performed by the attribution of different material properties to the Nucleus Pulposus and the Annulus Fibrosus. On Abaqus (Dassault Systèmes, USA), the model was assembled with the cervical spine ligaments. For the ligaments, two parameters were studied: insertion area and number of elements. The range of motion (RoM) of the model, measured on the superior surface of each IVD, was compared with experimental mobility data conducted in Panjabi et al. (2001). Based on the RoM, the influence of patient-specific geometrical features was evaluated.

For validation, the FE model was conditioned to a simultaneous 50N compressive load and 1Nm moment in different directions. The FE model agreed well in flexion where motion for each segment was 7.97° at C4-C5, 5.58° at C5-C6 and 2.84° at C6-C7, as well as in extension where motion for each segment was 6.80° at C4-C5, 4.65° at C5-C6 and 1.67° at C6-C7 which are within the experimental findings. For lateral bending, motion for segment at C4-C5 was 7.92° and at C5-C6 was 5.59° were in accord with the experimental data, however motion for segment at C6-C7 was 2.15° which meant an underprediction by 153%. For axial rotation, motion for segment at C4-C5 was 9.48° and at C5-C6 was 7.35° which overpredicted by 29% and 28%, respectively. For segment at C6-C7 motion was 3.60° which is within the experimental data.

Adjustments on the cervical spine components produced favorable outcomes, when compared to the experimental mobility data. The workflow followed in this study ensured uniform segmentation of cervical spine finite element models and latter adaptations of parameters. Compared to non-linear FE models, this methodology presented similar mobility results with a facilitated segmentation process and linear material properties.

Bone tissue engineering (BTE) scaffolds are used in different situations when there is a need to replace or regenerate bone tissue, such as fractures, tumor resections or joint replacement. Among biocompatible materials bio-Ceramic scaffolds represent an effective solution for BTE applications since they stimulate the bone growth thanks to their affinity with the native bone tissue. The main drawback of ceramic materials is represented by their intrinsic brittleness that makes the failure unpredictable and catastrophic. One of the key challenges in tissue engineering is to design scaffolds with appropriate mechanical properties and pore structures that can carry the load and in the meantime support cell growth and differentiation. To this purpose, we developed a numerical tool able to generate bone-like architectures with tuneable morphometric properties based on stochastic approaches. In particular, the Voronoi tessellation was used to create the three-dimensional scaffold mimicking the natural structure of bone tissue. A predefined number of seeds have been randomly distributed in a cubic domain and the Voronoi tessellation have been created. According to a statistical probability function based on the orientation of the edge of the cells, faces and edges have been removed in order to get a fully interconnected porosity. Finally, the resulting structure have been dilated to get the scaffold architecture. Taking inspiration from the trabecular bone tissue, the model was forced to have a preferential orientation of the rod-like trabeculae along a prescribed direction, no more than four trabeculae converging in one node and control the ratio between rod-like and plate like trabeculae.

The homogenized elastic properties of the scaffolds were characterized using finite element analysis on cubic domains by using the solver ParOSol. A damage-based iterative procedure has been used to assess the macroscopic strength of the scaffold.

The proposed model offers a versatile and effective method for designing ceramic scaffolds for tissue engineering applications. The ability to control the scaffold properties through the Voronoi tessellation allows for the development of customized and patient-specific scaffolds with enhanced mechanical properties, potentially leading to more effective tissue regeneration and repair. In particular, through this model we were able to control the amount of plate-like trabeculae and rod-like trabeculae, their orientation and detect their effect on the macroscopic strength of the scaffold.

In the field of bone tissue engineering (BTE), porous scaffolds gained attention due to their capability to promote tissue regeneration. Recently, triply periodic minimal surfaces (TPMS) scaffolding method became a promising candidate to aid bone tissue renewal. TPMS scaffolds allow tuning of several designing parameters, so the aim of this study was to understand the role of geometry, porosity degree and material non-linearity in the overall mechanical response of the scaffold. For this purpose, different building materials for TPMS scaffolds were simulated using the finite element (FE) method.

The TMPS scaffolds derived from a combination between geometry, Diamond (SD), Gyroid (SG) or Primitive (SP); and porosity, 60%, 70% or 80%. The nine different scaffolds (SD60, SD70, SD80, SG60, SG70, SG80, SP60, SP70 and SP80) were modelled on Abaqus (Dassault systems Simulia Corp., USA). The FE models were generated by the repetition of a cubic structural unit (1 mm side) on a 2x2x4 units configuration and modelled using two materials: VisiJet^® M3 Crystal (approximated as linear elastic, with E=1.46 GPa and ν=0.35), and a viscoelastic gelatin hydrogel (Kalyanam et al., 2009). On a first stage, the mechanical behaviour of the 16-units FE model was compared to a previous study that used FE models of 2-units (Castro et al., 2020). Preliminary results showed different mechanical behaviour on the maximum value and location of von-Mises stress, between 16- and 2- units. Thus, the 16-units scaffolds were considered for the following analysis. The mechanical loading was performed as a displacement-controlled compression along the z-direction of the FE model, applying a compression amplitude (measured along the z-direction) of up to 10% for the linear material and 2% for the hydrogel.

The scaffolds built with the linear elastic material showed a linear trend for each combination of geometry and porosity (R² coefficient from 0.99466 for SP80 to 0.99994 for SD80), in terms of equivalent strain/stress curves, revealing the absence of non-linear behaviour coming from the changes in scaffold geometry. In terms of equivalent stiffness, for each geometry, the Young’s modulus decreases with an increase in porosity, as expected. Furthermore, the SD and SP geometry corresponded to the highest and lowest Young’s Modulus, respectively, for each porosity level. The computational equivalent stiffness was compared with the stiffness values obtained from experimental tests (3D printed scaffolds in Crystal) and homogenisation method, showing considerable differences for the SD geometry and minor variations for SG and SP. However, the homogenisation method assumes an infinite periodicity of units, which probably means that, by increasing the number of units on the FE model, the outputs from these two methods could be closer.

The use of computational FE models allowed to study the mechanical behaviour under compression of TPMS scaffolds. The current outputs suggest that the non-linearities on the mechanical behaviour of TMPS scaffold are not introduced by the scaffold geometry, allowing for the exploration of non-linear building materials. Further work will include experimental tests of stress relaxation, to investigate the effect of geometry on the time-dependent behaviour of these scaffolds.

The application of computational methods to the design and testing process of Bone Tissue Engineering scaffolds is a valuable but complex approach, as the number of structures and materials that can be utilized to create these structures is constantly increasing and with it the time and costs required to accurately assess scaffold properties.

Among the different parameters that need to be defined, the macroscopic strength and the fracture mechanisms are of particular importance for ceramic-based scaffolds, as they are characterized by brittle-like failure and can undergo abrupt rupture processes. Proper assessment of scaffold strength and fracture patterns can aid the design and optimization process and therefore needs to be accurately implemented in the scaffold testing process.

Different crack propagation models exist nowadays, allowing for the computational evaluation of scaffold strength. Explicit methods can be characterized by stability issues, as convergence is conditionally achievable, and tends to limit its analyses to short transient processes. The phase-field approach, which is based on an implicit approach, models crack propagation as a gradient of element damage, allowing for the evaluation of brittle fracture as a quasi-static process. Furthermore, its versatility allows for the implementation of different phenomena, such as corrosion.

The staggered phase-field model created by Molnár et al. for modelling brittle failure was tested on hydroxyapatite-based scaffolds produced with the Robocasting technique. The results of this starting algorithm highlighted discrepancies in the failure mechanism when compared to finite element method results. While the latter displayed vertical cracks along the compression direction, suggesting that tensile stress is the dominating factor for fracture generation, the phase-field results were characterised by diagonally oriented cracks, highlighting a shear-dominated fracture process. This particular failure mode is owed to the specific function assumed for onset and progression of localised damage.

In this work, an alternative form of damage propagation is proposed with the purpose of better simulate the fracture mode observed in ceramic scaffolds subjected to compressive loads.

The aim of this alternative formulation is to account for a higher resistance of the material to shear stresses, and a much higher resistance to compression.

The final crack pattern formation more closely resembles the expected pattern and is therefore able to model the material’s behaviour more accurately, paving the way to the application of the phase-field model on bone tissue engineering scaffolds characterized by structural continuity between different layers along the loading direction, where compression stresses are transferred directly along said direction.

Skin is the body’s largest organ, composed of three layers that vary in properties across its structure and exhibit directional-dependent physical properties. Despite its exceptional healing capacity, the skin is very susceptible to full-thickness wounds that require specialized dressings to promote tissue regeneration. While there are various commercial skin grafting solutions available, they all have limitations that can result in non-native scar formation, for instance, highlighting the need for further investigation in the field. In this study, we investigate the use of fiber scaffolds with micro-scale auxetic geometries for skin tissue engineering applications. Auxetic structures are known to expand in multiple directions when stretched, enabling them to conform to different wound shapes and promote wound healing [1]. By altering the micro-geometry design, these structures can be tuned to provide tailored mechanical properties that match the unique needs of each patient’s wound, which may vary based on location of the wound or patient’s age or gender [2].

We have investigated two types of scaffold micro-scale designs, both comprising interwoven fiber networks to achieve auxetic behavior. These designs are categorized as re-entrant and chiral and exhibit distinct mechanical properties that are determined by the arrangement of interconnections between the fibers. Numerical finite element simulations of these fibrous structures were conducted to predict their mechanical behaviour prior to fabrication by means of a parametric study of each auxetic design. These simulations were generated, calculated and post-processed with a Python script-based software tool developed to automate this process. After numerical analysis, some designs were selected for fabrication and mechanical characterization with uniaxial and biaxial tensile tests to validate the numerical models. Polycaprolactone fiber scaffolds were fabricated using an in-house built melt electro-writing set-up. The multilayer structures were printed with fiber diameters of approximately 20 microns (with pore size around 300 microns) with the aim of developing microstructures that could provide adequate mechanical environment for cells.

Both numerical simulations and experimental tests have revealed that the mechanical behaviour of the printed scaffolds matches the mechanical properties [3] and behavior of skin (J-shaped stress-strain curve), as well as its auxetic behaviour [4]. Mimicking native tissue properties is essential for tissue engineering. Furthermore, comparative analysis of the results revealed differences in stiffness and auxeticity depending on the auxetic design, allowing us to tailor their properties via geometrical changes.

The combination of numerical simulations with a high-end fiber printing process and mechanical testing has enabled us to develop a platform for designing structured fiber scaffolds with tunable mechanical properties for skin tissue regeneration.

References

1. Kim et al, Materials (Basel), 14(22):6821, 2021.
2. Luebberding et al, Skin Res Technol, 20(2):127-35, 2014.
3. Joodaki et al, Proc Inst Mech Eng H, 232(4):323-343, 2018.
4. Dwivedi et al, R Soc Open Sci, 9(3):211301, 2022

Acknowledgements

The authors acknowledge the project LMP 176_21 funded by the Department of Science, University and Knowledge Society of the Government of Aragon.

The complex process of bone remodelling is strongly linked to the electric fields that naturally occur in the bone as a bioelectric tissue. Still, this process needs to be fully understood. Mechanical loading causes charge shifts in the bone, which are characterised by streaming potentials and the piezoelectric properties of the collagen matrix of the bone.

The biological structure of bone is multiscale. The entities on different spatial scales, such as osteoblasts/osteoclasts on the microscale, bone tissue on the mesoscale, and whole bone on the macroscale, show responses to endogenous electric fields in distinct time scales. Electrical stimulation aims to imitate these naturally occurring electrical signals to enhance bone regeneration. The design of innovative implants that employ electrical stimulation for accelerated bone regeneration is a fast-growing research area. In recent years, much research has been done to develop electrically stimulating implants for different applications, e.g., bone fractures, non-unions, or avascular necrosis. The stimulation parameters used in clinical and in vivo applications are mainly used empirically, and little is known about the actual field distribution inside the tissue. Numerical simulation can close this gap and help predict the tissue's electric field distribution, optimise stimulation parameters, and identify essential parameters, e.g. via Uncertainty Quantification.

On the macroscale, we will present different examples of finite element models of electrostimulating bone implants for bone defects, e.g., in the jaw region. More specifically, we will discuss practicable stimulation arrangements and electrode designs and optimise stimulation parameters (dimensions of the electrodes, stimulation voltage) with regard to the electric field strength as a target variable. Modelling parameters with an essential impact on the electric field distribution will be identified via sensitivity analysis. Furthermore, on the mesoscale, the electric fields experienced by the bone cells and their response affect bone remodelling. On this front, we present a mathematical model based on the mean-field analysis, which explains the effect of the electric field on bone cells in an in-vitro setup. In addition, on the microscale, each bone cell changes its shape and adhesion with the orientation of the electric field. The cell contractility affecting its shape and adhesion can be mathematically modelled as a multi-physics phenomenon, particularly a bio-chemo-electro-mechanical model.

With these phenomena that happen in distinct time scales across the spatial scales of the bone, it is necessary to implement a multiscale numerical simulation bridging these scales. Such a multiscale model helps provide a better understanding of the impact of electrical stimulation on bone regeneration. Further, such models can enhance the experimental characterisation of therapeutical electrical stimulation of bone. In this context, we present our approach to developing a multiscale framework for electrically active bone implants.

Bone can optimize its structure under the influence of mechanical loads. The macroscopic load acting on bone influences the bone cells, particularly osteocytes present in the lacunae canalicular network (LCN). These osteocytes respond to numerous physical signals, including substrate strain, interstitial fluid flow, and pore pressure. Cellular-level strains in vivo are considered too small to elicit a response from osteocytes. However, several studies suggest that load-induced interstitial fluid flow (IFF) in LCN may be a primary stimulus as it exerts shear and drags forces on osteocytes. Besides IFF, it has also been shown that pore pressure generated under physiological loading conditions is adequate to enable osteocytes to respond. Despite the importance of IFF and Pore pressure, relatively few studies explore their influence on bone adaptation. Motivated by the fact above, we investigated the role of interstitial fluid flow and fluid pore pressure on osteogenesis.

This work aims to predict new bone formation on both cortical bone surfaces under mechanical loading. We use dissipation energy density due to fluid flow and fluid pore pressure as a stimulus because, being a scalar quantity and having attributes similar to fluid motion, it is convenient to use. We hypothesize that the site-specific new bone distribution at the endocortical and periosteal surfaces is directly proportional to the square root of the summation of the dissipation energy density due to canalicular fluid flow and pore pressure. Accordingly, a poroelastic finite element model of a simplified geometric beam with a cross-section similar to the mid-cross-section of 16 weeks old C57BL/6J was subjected to the mid-strain axial loading protocol Barmen et al. [1]. The fluid velocity and pore pressure estimated from the above analysis is used to calculate the dissipation energy density. This model is coupled with a novel mathematical formulation, considering the summation of the dissipation energy density of poroelastic flow and pore pressure as a stimulus to predict cortical bone adaptation.

The bone formation rates (BFR) calculated by the model are 0.9986 µm³/µm²/sec and 0.6348 µm³/µm²/sec, respectively, at the periosteal and endocortical surfaces, which are not significantly different from the corresponding experimental BFR of 0.9257 ± 1.5855 µm³/µm²/sec (p-value = 0.9817, t-test) and 0.5964 ± 1.4285 µm³/µm²/sec (p-value = 0.9791, t-test). Additionally, the statistical significance of site-specific MAR was measured using Watson’s U² test. The computed MAR at both the cortical surfaces was not significantly different from the experimental MAR at the periosteal surface (p-value = 0.3417, Watson’s U² test) and at the endocortical surface (p-value = 0.94, Watson’s U² test).

Reference:

[1] A. G. Berman, C. A. Clauser, C. Wunderlin, M. A. Hammond, and J. M. Wallace, 2015, PLoS One, 10 (6), 1-16.

Fracture healing is a complex process that involves numerous biological and mechanical factors. Understanding this process is crucial for healthcare professionals involved in the treatment and management of fractures, but also, for students pursuing careers in fields such as medicine, physical therapy, and biomechanics. Learning about fracture healing can provide them with a foundation for understanding musculoskeletal injuries and developing effective interventions and implant designs. In addition, gaining knowledge of fracture healing can inspire future research and innovation in the field, leading to improved treatments and outcomes for patients. However, with advancements in technology and increasing accessibility of simulation models, there is potential for a wider range of individuals to gain hands-on experience through simulation and digital learning methods. But up to date, simulation of fracture healing is limited to research institutes and highly trained professionals.

On the basis of the Ulm fracture healing model (Shefelbine et al. 2005, Simon et al. 2011), we developed a software platform (OSORA Medical Fracture Analytics), which is able to be used for educational purpose as well as for research questions.

We used the software in a course on fracture healing, where computational science and engineering students were exposed to different questions as the influence on the healing performance of:

• Fracture gap size
• Fracture angle
• Osteosynthesis stiffness

Into the web-based tool, we implemented a wizard to generate simulation input by students in an intuitive manner. Adjustable parameters where: AO-fracture class, geometrical bone and fracture dimensions, osteosynthesis material and loading conditions. Simulations were then automatically started and results available online.

The learning objectives for the students were the research of suitable literature data on experiments and studies to corroborate the simulation results, simulation performance and result interpretation.

Deriving qualitative and quantitative results from the web-based tool allowed a detailed comparison with literature data on those influences for the students.

The influence of fracture gap size was compared with in vivo data from Claes et al. 1997, Meeson et al. 2019 and Markel/Bogdanske 1994. Trends in the influence of the fracture angle where compared with Park et al. 1998 and Yamagishi et al. 1955. Literature data of simulations from Steiner et al. 2014 were used to explain and compare the influences of intramedullary nail stiffness behavior on the healing outcome.

Quantities as the interfragmentary motion, bone formation tissue volume over time and longitudinal fracture stiffness where analyzed and compared to the experimental findings.

With this web-app the students were able to in deep discuss on healing patterns and quantitative comparisons with literature data. The platform offers a flexible and engaging learning experience that was integrated into traditional classroom instruction, providing an efficient and effective means of enhancing biomechanical and medical education with computer simulation.

Introduction: While Wolff’s law established bone’s adaptation to mechanical environment, Frost’s mechanostat theory predicted a mechanical threshold (such that strain, strain energy density etc.) for such adaptation [1,2]. Such threshold (for example, strain) was, however, found to vary with forcing frequency [3]. Moreover, strain or strain energy density alone does not predict new bone formation [4]. Evidence of involvement of interstitial fluid flow in bone adaptation has been increasing in the literature [5]. Accordingly, bone adaptation models based on poroelastic dissipation energy density have been developed by various researchers [6-8]. The present work advances the poroelastic model by establishing the strain threshold and bone formation rate (BFR) as a function of forcing frequency.

Methods: Biot’s poroelastic model has been studied in frequency domain, and dissipation energy density has been calculated for a sinusoidal forcing frequency using standard methods [9]. Data from the literature have been used to establish relation among strain threshold, BFR and forcing frequency [3].

Results: BFR was found to be proportional to “square-root of dissipation energy density” minus a corresponding threshold value (of the “square-root of dissipation energy density”). The expression is similar to that obtained by the authors in previous study based on different experimental data, which were for non-sinusoidal waveform of loading on mouse tibia [8]. Dissipation energy threshold was found to be constant that does not depend on forcing frequency. On the other hand, the corresponding strain threshold was found to be a function of frequency. There is a frequency predicted for which the strain threshold is the minimum and that frequency is around 20 Hz, which agrees with the literature [10]. The strain threshold needed for new bone formation monotonously increases beyond this forcing frequency.

Conclusions: This work establishes that the rate of change of BFR with respect to square root of dissipation energy density is a constant, regardless of forcing frequency. The strain threshold has also been expressed in terms of frequency, which can easily determine the simple sinusoidal loading regimen (peak strain and frequency) to maintain a bone mass and mitigate bone loss arising due to bed-rest, space-flight, muscle paralysis etc.

References:
[1] Wolff, 1870, Arch. Für Pathol. Anat. Physiol. Für Klin. Med. 50, 389–450.
[2] Frost, 1987, Anat. Rec. 219, 1–9.
[3] Hsieh & Turner, 2001, Journal of Bone and Mineral Research, 16(50), 918-924.
[4] Tiwari & Prasad, 2017, Biomech. Model. Mechanobiol. 16, 395–410.
[5] van Tol, Schemenz, Wagermaier, Roschger, Razi, Vitienes, Fratzl, Willie, & Weinkamer, 2020, Proc Natl Acad Sci, 117(51):32251-32259.
[6] Kumar, Jasiuk, & Dantzig, 2011, Journal of Mechanics of Materials and Structures, 6(1-4), 303-319.
[7] Pereira, Javaheri, Pitsillides, & Shefelbine, 2015, J. R. Soc. Interface, 12, 20150590.
[8] Singh, Singh, & Prasad, 2021, Proceedings of the ASME IMECE2021, 71220.
[9] Biot, 1941, Journal of Applied Physics, 12, 155-164.
[10] Rubin, & McLeod, 1994, Clin Orthop, 298, 165–174.