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.

Different computational-based approaches have been extensively used for simulating the morphogenesis of multicellular systems. On the one hand, continuum-based models allow the simulation of large cell populations at the macroscopic level. This is the case for reaction-diffusion systems based on partial differential equations (PDEs) or positional information (PI). On the other hand, discrete approaches with agent-based models considering cells as autonomous entities that interact among themselves and with the microenvironment. These models have been widely employed, for instance, to simulate tumor growth in vitro [1], to study the role of extracellular matrix density in cell migration within solid tumor spheroids [2] or lumen-based organoids [3], collective cell migration [4].

In this work, we present the results of different in-vitro experiments corresponding to the morphogenesis of different human-cell-based tumor organoids, from neuroblastoma, pancreas and lung. Normally, all these experiments follow a similar protocol [5], where individual tumor cells are cultured in 3D microfluidic devices upon collagen-type I hydrogels. With the time these tumor cells are able to proliferate, self-organizing according to the architecture of the hydrogel, creating different structures. The evolution of size and shape of these tumor structures are quantified by means of microscopy-based images. Thus, this kind of organ-on-a-chip experiments constitutes a novel modelling strategy of in vitro multicellular human systems that in combination with numerical simulations provide a relevant tool for research in mechanobiology.

Therefore, in this presentation I will show our engineering-based strategy in which we integrate computational models and in-vitro experiments in order to tackle of how matrix mechanics is regulating the size and organization of tumour organoids. Different computer-based modelling strategies will be presented to simulate the morphogenesis of these in-vitro tumor organoids, clearly identifying the challenges that we face in these numerical simulations, mainly associated to the variability observed in the experiments and the difficulties associated to the parameters calibration process.

References:

[1] Hervas-Raluy, S., Wirthl, B., ... & Wall, W. A. (2023). Tumour growth: An approach to calibrate parameters of a multiphase porous media model based on in vitro observations of Neuroblastoma spheroid growth in a hydrogel microenvironment. Computers in Biology and Medicine.

[2] Gonçalves, I. G., & Garcia-Aznar, J. M. (2021). Extracellular matrix density regulates the formation of tumour spheroids through cell migration. PLoS computational biology, 17(2), e1008764.

[3] Camacho-Gómez, D., García-Aznar, J. M., & Gómez-Benito, M. J. (2022). A 3D multi-agent-based model for lumen morphogenesis: the role of the biophysical properties of the extracellular matrix. Engineering with Computers, 38(5), 4135-4149.

[4] González-Valverde, I., & García-Aznar, J. M. (2018). Mechanical modeling of collective cell migration: An agent-based and continuum material approach. Computer Methods in Applied Mechanics and Engineering, 337, 246-262.

[5] Plou, J., Juste-Lanas, Y., Olivares, V., Del Amo, C., Borau, C., & García-Aznar, J. M. (2018). From individual to collective 3D cancer dissemination: roles of collagen concentration and TGF-β. Scientific reports, 8(1), 12723.

Acknowledgments: This work was supported by the European Research Council (ICoMICS Adv grant: 101018587 and the Spanish Ministry of Science and Innovation (PID2021-122409OB-C21).

Physical forces are critical for the development, maintenance, and adaptation of various tissues in the body including bone. Bone adapts to mechanical loads through actions of bone cells: osteoclasts that specialize in bone resorption; osteoblasts that produce bone; and osteocytes embedded in bone matrix that sense mechanical forces and coordinate bone adaptation. ATP is a ubiquitous intracellular molecule critical for cellular bioenergetics, which is also one of the first signals released from bone cells that experience mechanical stimulation. Extracellular ATP is degraded by ecto-ATPases to form ADP and eventually adenosine. We developed a theoretical model describing mechanically induced ATP and ADP release, followed by their extracellular diffusion and degradation. The model was validated using experimental data for calcium signaling, ATP and vesicular release upon mechanical stimulation of a single osteoblast, and then scaled to a tissue-level injury. The route and amount of ATP release depended on mechanical forces, ranging from vesicular release of small ATP boluses upon membrane deformation, to leakage of ATP through resealable plasma membrane tears, to spillage of cellular content due to destructive forces. The total amount of ATP released indicated the degree of injury and determined the maximal distance from the injury at which responses were stimulated. The peak ATP concentration discriminated between injuries that released similar amounts of ATP due to differences in cell repair, and determined signal propagation velocity. ADP-mediated signaling was relevant for larger tissue-level injuries, conveying information about the proximity to the injury site and its geometry. After their release from an injury site, ATP, ADP, and adenosine signal through purinergic receptors, including seven P2X ATP-gated cation channels, seven G-protein coupled P2Y receptors responsive to ATP and ADP, and four P1 receptors stimulated by adenosine. All bone cells express multiple P2 receptors. We used published studies of overexpressed individual P2 receptors to fit their dependencies on ATP concentration using the Hill equation, and experimentally examined the concentration dependences of P2-mediated calcium responses to ATP and ADP in osteoblasts/osteocytes and osteoclasts. Dependence of the calcium response on ATP concentration exhibited a complex pattern that was not explained by the simple addition of individual receptor responses. To deconvolute the contribution of individual receptors to the ATP dose-response curve, we constructed a mathematical model that combined the actions of the P2Y2 receptor that has the highest affinity to ATP and induces inositol trisphosphate (IP3) mediated calcium release, and the P2X7 receptor that has the lowest ATP affinity and allows calcium influx through the pores. The model predicted interactions between these P2 receptors at the levels of receptor interactions and cellular signaling, which were validated using CRISPR/Cas9-generated P2Y2 and P2Y7 knockout cells. Our studies demonstrate how spatiotemporal signals provided by extracellular ATP and ADP encode the information regarding severity, scale and proximity of the mechanical stimulus, allowing for a well-choreographed tissue responses to physical forces.

1. Introduction

Mechanobiology is an area of intense research focus nowadays. In cellular length, it tries to describe how mechanical forces can modulate the structure and function of cells. These forces are transmitted to the nucleus, through direct physical connections to the cytoskeleton via LINC [1], and deform the nucleus, which in turn triggers different biochemical pathways. Therefore, nuclear deformations are critical structural and biophysical parameters in mechanobiology [1].

Recent studies have highlighted the variations of cellular responses with the surrounding stiffness [2]. Inspired by this, we decided to use the finite element method to investigate variations of nuclear deformations with different mechanical environments. For this aim, we modeled the realistic geometry of a single cell in Abaqus^® to measure to deformation applied to the nucleus after the contraction of the acto-myosin units.

2. Materials and Methods

We created a 3-D model of a cell with a volume of 537 nm³ including the membrane, nucleus, nuclear cap, dorsal and ventral f-actins, LINCs, and cytoplasm. We assumed that the dorsal and ventral f-actins are tied to the cell membrane in both extremities. Therefore, fibers contraction would deform the cell membrane and the nucleus as well due to, first, the contact forces between fibers and the nuclear cap, and second, the forces coming from the LINCs. Moreover, to take into account the interaction of the call with its surrounding environment, we attached it to the top surface of a substrate to simulate a traction force microscopy test. As shown in figure 1, the modeled cell has been tied to the substrate in 3 focal adhesion sites [1].

3. Results

We measured the deformation applied to the nucleus for a substrate with 3 different Young’s moduli: 2, 4, and 8 kPa. The average values of the maximum principal strain in the nucleus are 5.3%, 5.6%, and 6.0%, respectively. Figure 2 also shows the distributions of the maximum principal strain for these 3 cases. These results suggest that there is a direct relationship between strain values in the nucleus and the stiffness of the surrounding environment.

4. Discussion and Conclusions

We developed the first finite-element model of a realistic and complex single cell to quantify deformations of the nucleus after the contraction of the acto-myosin units. We showed that finite element modeling can also be regarded as a strong tool to measure the interactions between different cellular components. The simplicity of the proposed model is a promising asset for further quantifications of intracellular mechanics in 2-D and 3-D environments.

5. References

1. Ghosh S et al., Cell reports, 27(5): 1607-1620 (2019).

2. Petit C et al., Biomechanics and Modeling in Mechanobiology, 20(2):717-731 (2021).

The cortex is a thin layer located beneath the membrane of animal cells. It governs important cell mechanical properties, including cell shape. It is composed of a dense network of actin filaments and actin-binding proteins, in particular myosin motors. Myosins generate active stresses in the actin network and make it behave as an active material. The spatiotemporal organization of the cortex is tightly coupled to that of signaling molecules, in particular Rho GTPases. Indeed, the interaction between actin, myosin and Rho GTPases can give rise to self-organized patterns. Thus far, studies have focused on spatially extended patterns, like actin polymerization waves. However, spatially localized patterns are also observed in the cortex, whereby isolated spots enriched in actin, myosin and signaling molecules play a role in crucial processes, including cancer invasion. What is the origin of such spatially localized patterns? We use a simple physical theory to show that they can emerge from the coupling between active cortical mechanics and signaling reactions, through an instability called "slanted snaking". Beyond single cells, our theory could also explain the origin of spatially localized, mechanochemical cues in tissue morphogenesis.

Cytoskeletal networks play a key role in multiple mechanical and dynamical processes in cells. Recently, a continuum theory has been developed [1], allowing for the prediction of the material properties of highly crosslinked cytoskeletal networks from a phenomenological modelling of the microscopic interactions between the cytoskeleton filaments. We extend this theoretical framework to account for external forces, allowing us to explore how the properties of cytoskeletal networks are affected by the presence of various interfaces. This extended theory allows to study the interplay between cytoskeletal networks and the organelles that it interacts with,such as the centrosomes in spindles, vesicles embedded in intracellular actin networks, or even cortex-membrane interactions at the cell periphery.
[1] S. Fürthauer, D. J. Needleman, M. J. Shelley, NJP 23, 013012 (2021).

Eukaryotic cells exhibit a complicated rheology in response to mechanical stimuli, including an elastic response due to the cell cytoskeleton (a network of crosslinked filamentous proteins) and energy dissipation resulting from transport of cytosol through this network as well as transient crosslink dynamics [1]. Existing models of cytoskeletal mechanics fall into two categories: discrete (microscale) models enable inclusion of detailed biophysics but are typically computationally challenging due to the large number of discrete elements and their interactions, while continuum (macroscale) models are easier to solve but the manner in which microscale parameters and processes manifest themselves at the macroscale is usually unclear. Mathematical modelling efforts involving more rational and rigorous mathematical methods (such as discrete-to-continuum upscaling or homogenization) to systematically bridge between these two approaches are still largely missing. Our aim is to develop a multiscale framework providing such a bridge.

We introduce a discrete mathematical model for the mechanical behaviour of the eukaryotic cell cytoskeleton during intracellular transport. The model involves an initially regular array of pre-stretched protein filaments (e.g. actin, vimentin) which exhibit resistance to enthalpic stretching, joined at crosslinks to form a planar network. To mimic inertialess motion of a small object placed in the domain, we impose a quasi-static displacement of a set of crosslinks in the centre of the domain, solve for the remaining nodes through a local force balance, and calculate the net force on the object. Assuming that the inter-crosslink distance is much shorter than the lengthscale of the cell, we rationally upscale the force balance equations using the discrete-to-continuum method based on Taylor expansions to form a continuum system of governing equations, inferring the corresponding macroscopic stress tensor and strain energy.

We solve these discrete and continuum models numerically and infer force-displacement curves, which show good quantitative agreement across the parameter space. The force-displacement curves have weak dependence on the pulling angle (with respect to the initial filament orientation). The net force acting on the object increases with increasing pre-stress and larger objects. Furthermore, we linearize the continuum model to construct analytical approximations for the stress and strain fields in the neighbourhood of the moving object, and explicitly compute the net force required to generate a given deformation as a function of key parameters, such as the object and mesh sizes, the pulling angle and the network pre-stress.

Our mathematical formulation allows us to make explicit predictions of the force-displacement curve in optical-tweezers experiments [1], analytically characterizing the (linearized) rheology of the cytoskeleton. Future work will also incorporate nonlinear effects in polymer elasticity and dynamic aspects of cell rheology (poro-visco-elasticity). We also plan to apply similar ideas to other biopolymer networks, namely collagen fiber networks forming extra-cellular matrices.

Acknowledgements:

This research was funded by EPSRC grant EP/S030875/1. We thank Profs. Ming Guo and
Roger Kamm (MIT) for valuable discussions.

References:

[1] Hu, J. et al. “Size-and speed-dependent mechanical behavior in living mammalian
cytoplasm.” Proceedings of the National Academy of Sciences (2017): 9529-9534.

Abstract:

Cardiovascular diseases are the major cause of death around the world. It is estimated that 20% of the population is affected by elevated arterial wall stiffness which increases the risk of aneurysms rupture. Also, in the case of narrowed vessels due to arterial wall remodelling the blood flow is disturbed which may lead to an increased hemodynamic pressure gradient and increased cardiac load. Therefore, the parameters describing the wall stiffness and pressure gradient are used as indicators helping with diagnosis of disease severity [1]. Computational models of hemodynamics such as lumped parameter models (referred to as 0D), 1D models and 3D Fluid-Structure Interaction models (FSI) are tools for studying the cardiovascular system and may be further applied to non-invasive diagnostics and disease research [2, 3].

The parameter estimation problem may be solved by variational or sequential approach. The sequential approach is iterative in nature and requires many simulations of the forward problem which may be prohibitive in the case of FSI simulations. With a sequential approach based on the Kalman Filter, the model prediction is improved at every time step by measuring the discrepancy between model output and measurements. It was shown that the total computational time for the sequential approach is of the same order of magnitude as the CPU time needed for one forward simulation [4].

For FSI models input parameters like wall stiffness and boundary conditions strongly affect the solution. Therefore, the challenge in constructing a model is the determination of parameters which make the simulation results agree with clinical data. The Kalman filter and its various modifications are used as methods to estimate unknown parameters that may not be directly observable [2].

The goal of this work is to explore the applications of Kalman filters to estimate the unknown parameters and boundary conditions for hemodynamic models and apply it to 1-way and 2-way FSI model of the carotid artery to estimate the Young’s modulus.

Acknowledgments:

The research leading to these results is funded by the Norwegian Financial Mechanism 2014-2021 operated by the National Science Center, PL (NCN) within GRIEG programme under grant UMO 2019/34/H/ST8/00624, project non-invasivE iN-vivo assessmenT Human aRtery wALls (ENTHRAL, www.enthral.pl)

[1] D. Nolte, C. Bertoglio, Inverse problems in blood flow modeling: A review, International Journal for Numerical Methods in Biomedical Engineering 38 (8) (2022) e3613.

[2] C. J. Arthurs, N. Xiao, P. Moireau, T. Schaeffter, C. A. Figueroa, A flexible framework for sequential estimation of model parameters in computational hemodynamics, Advanced modeling and simulation in engineering sciences 7 (1) (2020) 1–37.

[3] A. Quarteroni, A. Veneziani, C. Vergara, Geometric multiscale modeling of the cardiovascular system, between theory and practice, Computer Methods in Applied Mechanics and Engineering 302 (2016) 193–252.

[4] Bertoglio, Cristóbal, Philippe Moireau, and Jean‐Frederic Gerbeau. "Sequential parameter estimation for fluid–structure problems: application to hemodynamics." International Journal for Numerical Methods in Biomedical Engineering 28.4 (2012): 434-455.

Aortic root connects the left ventricle to the ascending aorta and houses the aortic valve (AV) ensuring one-direction flow of blood during systole. The AV is normally composed of three leaflets, known as tricuspid aortic valve (TAV), but 1-2% of the population is born with only two leaflets, known as bicuspid aortic valve (BAV). The patients with BAV are considered at high risk of developing aneurysms and eventually dissection. The biomechanics of aortic root tissues are hypothesized to play an important role in the disease development. In this study, we use in-vivo echocardiographic images from TAV and BAV patients to analyze the differences in the biomechanics of aortic root tissues.

3D transesophageal echocardiographic (TEE) images of the aortic root were retrospectively acquired from 16 patients with the approval of the Institutional Review Board at the University of Pennsylvania. The images were segmented, registered, and converted into a medial model as presented in a previous study. The medial models were remeshed with an quadrilateral elements. Two methods were used for the biomechanical analysis: 1) patient-specific 3D inverse finite element (FE) modeling, and 2) population-level Bayesian inference based on radius variations.

The two approaches provided distinct advantages. The first, patient-specific approach preserves geometric details, but the effect of diastolic pressure and opening angle could not be accounted form. The second, Bayesian approach allowed us to calculate the population-level differences between TAVs and BAVs, but it discarded part of the information available from the images. The biomechanical differences we found in this work indicate that the aortic root tissue in BAV patients experience different intramural stresses that might be linked to the higher risk of aneurysm development. Future work will include implementation of growth and remodeling framework to further establish this link.

The primary function of arteries is as conduits to allow the heart to efficiently deliver blood throughout the entire body. The stiffness of arteries is a key functional parameter that can alter this efficiency, and increased arterial stiffness is a reliable predictor of cardiovascular risk [1]. However, directly determining the stiffness of the arterial wall is essentially impossible in vivo, and proxy measurements such as pulse wave velocity and total arterial compliance are the most feasible clinical measures of arterial stiffness. These, however, only reflect the average stiffness of a region of the arterial network and do not provide information about the local stiffness of arteries. As arterial stiffness is determined by changes in the tissue composition at local levels and diseases like atherosclerosis, aneurysms and dissections occur in relatively localized regions of the arteries, methods to provide accurate information about the local material properties are desirable to enable further research and novel clinical approaches.

We present our implementation of an inverse solver for estimation of local arterial stiffness with a 1D fluid-structure-interaction model of the artery. The model consists of a axisymmetric domain representing a human common carotid artery. The fluid is modeled as a Newtonian fluid with an assumed parabolic flow profile throughout the domain. We investigate two arterial wall models based on the theory of linear elasticity. The first derives from the application of Laplace’s law and the simplifying assumptions of a thin wall and is one of the most widely applied models for pulse wave propagation models of the arteries [2]. The second model is a novel implementation based on thick-walled cylinder theory. We implemented a least-squares procedure to estimate the Young’s modulus parameter in both models, and then evaluated this on experimental data collected from a laboratory phantom. Two sets of boundary conditions were compared. First, direct experimental data of measuremed inlet flow rate and outlet pressure were used. Second, a more general approach of applying a parameterized common carotid inflow and Windkessel outlet was applied. The Young’s modulus estimated with the thin walled approach is in general smaller than that from the thick walled approach. The different boundary conditions produce some what different time courses of pressure and flow as well as a difference in estimate Young’s modulus. Further work will compare the estimated Young’s modulus with stiffness determined through direct tension testing. Additionally, the method will be applied to estimate local arterial stiffness of the common carotid artery based pressure measured by applanation tonometry and flow and geometry by ultrasound imaging.

Acknowledgments

The research leading to these results is funded by the Norwegian Financial Mechanism 2014-2021 operated by the National Science Center, PL (NCN) within GRIEG programme under grant UMO-2019/34/H/ST8/00624, project non-invasivE iN-vivo assessmenT Human aRtery wALls (ENTHRAL, www.enthral.pl)

References

[1] Laurent et al., Eur Heart J; 27(21), 2588-2605 (2006).

[2] Boileau et al., Intl J Num Meth Biomed Eng; 31(10), (2015)

Arterial stiffness is an established biomarker of cardiovascular health [1]. By combining non-invasive measurements and computational models, arterial stiffness can be inferred through solving an optimization problem. However, the non-invasive measurements are hampered by measurement errors, and some parameter values in the optimization problem must be assumed which introduces additional uncertainties. To apply a novel computational model in clinical diagnostics, uncertainty quantification needs to be performed [2]. Model parameters which lead to a large variation in the model prediction can be identified through a subsequent sensitivity analysis. A large number of model evaluations are needed to estimate these sensitivity indices, thus, limiting the application of such analysis to computationally expensive models and allowing relatively few uncertain inputs. Using the Multifidelity Monte Carlo Method (MFMC) [3, 4], we estimate Sobol main and total effect sensitivity indices of a common carotid artery (CCA) 3D-fluid-structure interaction (FSI) model by offsetting the computational burden to computationally affordable 1D- and 0D- models. Computational resources are thus distributed over the three levels of fidelity such that a few 3D-FSI model evaluations ensure accuracy of the sensitivity indices while the lower fidelity models are leveraged to reduce the computational costs of the sensitivity analysis.

We will fist consider the situation where the CCA is modeled as an idealized, straight tube. The same geometric and material parameters as well as boundary conditions are applied to all models. At the inlet, a parabolic physiological flow rate and wave form is prescribed and at the outlet, a three-element Windkessel model mimics the downstream vasculature. In the 1D-model, the artery consists of nodes along a straight line, and in the 0D-model, the artery is represented with a resistor and a capacitor as an electrical analog. The arterial wall is modeled as a linear elastic material and blood is assumed to be a Newtonian fluid. Uncertain model parameters are the vessel diameter, arterial wall thickness, and material parameters for the arterial wall. The uncertainty and sensitivities of the pulse pressure, average pressure, and the diameter change are assessed. We will present the method and implementation we have developed for the CCA model and preliminary results comparing the sensitivity indices estimated through the MFMC approach with the ones estimated with the same computational budget through Monte Carlo simulation of the 3D-FSI model.

Time permitting, results will be shown for hyperelastic material models and patient-specific anatomies.

[1] Laurent et al., Eur Heart J; 27(21), 2588-2605 (2006).
[2] FDA, Assessing the Credibility of Computational Modeling and Simulation in Medical Device Submissions (2021).
[3] Qian et al. , J Uncertainty Quantification; 6(2), 638-706 (2018).
[4] Gorodetsky et al., J Comput. Phys. 408, 109257 (2020).

Acknowledgments:
The research leading to these results is funded by the Norwegian Financial Mechanism 2014-2021 operated by the National Science Center, PL (NCN) within GRIEG programme under grant UMO-2019/34/H/ST8/00624, project non-invasivE iN-vivo assessmenT Human aRtery wALls (ENTHRAL, www.enthral.pl).

Mechanical ventilation is a life-saving therapeutic tool for treating patients with impaired pulmonary function. Despite the clear benefits of mechanical ventilation, it can cause irreversible damage to the lung tissue. The main obstacles for more protective and individualized ventilation strategies are the still insufficient knowledge and understanding of complex lung mechanics in both healthy and diseased states, which is mainly due to the limited ability for in vivo measurement. To shed light into this issue, great efforts have been made in the past in the development of computational lung models. Even most of advanced work in this context only focus on investigating the effect of ventilation on tissue strains and stresses, while coupling to the pulmonary circulation is mostly neglected so far. This is despite the fact that the main function of the lungs, namely gas exchange, takes place through a dense network of pulmonary blood vessels in the alveolar walls. Hence, the coupling between the respiratory system and pulmonary circulation is crucial for getting more insight into the main purpose of ventilation: adequate oxygen supply and carbon dioxide release, while keeping the tissue at a healthy state.
In this contribution we therefore present a physics-based, coupled porous media approach to computationally model air flow, blood flow, and gas exchange in the human lungs. Motivated by the structure of the lungs, larger airways and blood vessels are modeled as discrete zero-dimensional (0D) networks that are embedded into a three-dimensional (3D), three-phase (air, blood and tissue) porous medium, representing the smaller airways, smaller blood vessels and lung tissue in a homogenized manner. Further, the respiratory gases, oxygen and carbon dioxide, are modeled as chemical subcomponents of air and blood with a suitable exchange model in the porous domain. To connect the homogenized and the discrete representations of airways and blood vessels, respectively, a 0D-3D coupling method is used, which allows a non-matching spatial discretization of both domains. The method couples fluid flow and species transport in these phases via an outflow condition from the tips of the discrete networks into the 3D porous medium and vice versa.
Such a comprehensive coupled approach allows us to study the complex interplay of tissue deformation and perfusion, and its effects on oxygenation and carbon dioxide release. Further, the underlying multiphase porous media model can easily be extended to include additional phases, so that pathological conditions such as water accumulation in pulmonary edema can be studied in future stages of the model. We consider our model to be a promising base for investigating clinically relevant questions, which might contribute to an improved treatment in respiratory care.

Mechanosensing of cells to the surrounding material is crucial for their physiological and pathological processes. The emergent dynamics of cells arise from a variety of interactions between cells and their local environment. However, materials design to guide cellular responses is largely ad hoc due to the lack of comprehensive modelling techniques for quantitative understanding. In this paper, we propose a computational model, that couples cell dynamics and cell-materials interactions, to study the mechanosensing of fibroblast cells seeded on different hydrogels.

In recent years, micromechanism-based theoretical modelling has been proposed to capture the essential biophysical characteristics, such as the generation of contractile forces in the cytoskeleton (CSK) and cell–substrate interaction. The contractile force of a cell is primarily generated by the intracellular stress fibers (SFs). SFs are formed via phosphorylation of myosin and polymerization of actin filaments. Actin filaments in cytoplasm are connected by α-actinin proteins to form actin bundles, which are crosslinked by myosin II proteins. The actin-myosin structure generates contraction forces in CSK by the ‘walking’ of myosin II along the actin filaments. Deshpande et al. (2007) proposed a bio-chemo-mechanical model to capture these key mechanisms.

In our continuum model, we consider the following important features of SF formation: (1) an activation signal is essential to trigger the formation of actin-myosin contractile units, (2) SF association rate is dependent on the activation signal and the dissociation rate is dependent on the contractile force, and (3) the dynamic contraction of SF is similar to the well-established muscle contraction model, which can be modelled by a modified Hill model.

For the cell-substrate interactions, we have considered both specific interactions (i.e., binding of membrane molecules to substrate ligands) and nonspecific interactions (e.g., van der Waals, electrostatic forces, hydrogen bonding and steric repulsion) (Bell et al., 1984). The assembly of active/specific focal adhesion (FA) is represented by the aggregation of integrins on the membrane and binding to substrate ligands. A thermodynamic model to represent the chemical equilibrium between these two integrin states (Deshpande et al., 2008; McEvoy et al., 2017), has been adapted in this work.

This coupled model allows us to predict the coupled effects of substrate stiffness and thickness on stress fiber formation and disassociation, and affinity integrin density. We also examine the effect of substrate on the cell-cell communications of fibroblast cells.

Our modeling results have revealed that a cell can sense its neighboring cell by deforming the underlying substrate. Our simulations also provide physical insights in the enhanced mechanosensing capacity of collective cells. The present modelling framework is not only important for profound understanding of cell mechanosensing, but also has the potential to guide the rationale design of biomaterials for tissue engineering and wound healing.

The physiological behaviour of the cardiovascular system is highly affected by the mechanical response of arterial segments, that is in turn dependent from both tissue histological architecture and the contractile tone of smooth muscle cells. The former depends mainly on the different amount and arrangement of constituents (mainly, elastin and collagen fibers), while the latter on chemical drivers of vasoactivity, such as nitric oxide (NO) and reactive oxygen species (e.g., ROS and PN). It is also noteworthy that the problem is highly multiscale since global hemodynamic conditions (e.g., heart rates, resistance of downstream vasculature) highly affect local flow conditions, and hence the local pressure field and the internal stresses affecting biochemical pathways governing vascular tone. Detailed high dimensional models (2D or 3D) can generally be used to simulate local hemodynamics of specific arterial sites, while the whole arterial tree is generally described through low dimensional descriptions (i.e., lumped 1D approaches).

This work presents a comprehensive multi-scale and multi-field computational framework that accounts for: i) a lumped 1D description of the macroscale arterial tree; ii) a continuum 3D model at the microscale of the local chemo-mechano-biological response of arterial tissues (accounting for passive and active tissue behavior); iii) biochemical-dependent vasoconstriction and vasodilation (the NO-ROS-PN biochemical chain), and biochemical-dependent tissue remodeling (the GFs- MMPs biochemical chain). Simulations from 3D chemo-mechano-biological models drive how parameters of the lumped description vary as function of segment dilation, as well as tissue histology and vasoconstriction.

The applicative case study investigates the relationship between arterial vasodilation and vasoconstriction with physical exercise. The obtained numerical results are consistent with available experimental data for normal and spontaneously hypertensive phenomena.

Growing trees respond to mechanical disturbances, resulting, e.g., from environmental forces or gravity, by forming so-called reaction wood. The latter is called compression wood in gymnosperms and tension wood in angiosperms. It enables the tree to control its posture by reinforcing and reorienting the axes of stems and branches, which is a key prerequisite for reaching large heights. The movement is due to asymmetric cambial activity resulting in eccentric growth and varying growth strains. While the underlying biological mechanisms of growth strain generation are not yet fully understood, various hypotheses correlating the induced macroscopic movement with the difference in cell wall structure of reaction and non-reaction wood have been proposed. On that basis, a homogenization procedure was developed for upscaling and evaluating the macroscopic effect of growth strains implemented at the cell wall level. In particular, we aim at unraveling the effect of the characteristic composition and microstructure of tension wood containing G-layers on the elastic properties and macroscopic growth strains by employing multiscale homogenization modeling techniques, based on the concept of continuum micromechanics. A three-step homogenization scheme is employed to estimate the elastic properties of the cell wall layers based on the volume fractions of elementary components. The structural organization at the cell wall level is represented by multilayered cylindrical inclusions exhibiting transversely isotropic material behavior taking into account the previously evaluated layer properties depending on the composition of their elementary constituents, such as cellulose, hemicellulose, lignin and water. We derive the stress and strain fields corresponding to four different loading conditions, which allows for constructing the complete stiffness tensor of tension wood and upscaling the growth strains induced within the G-layer. Coupling the micromechanics model with non-linear beam mechanics, applied to a growing branch inclined with respect to the vector of gravity, allows to simulate the reorientation process induced by growth strains at the cell wall level. In combination with experimental data of the branch shape evolution found in literature of specific species, growth-related parameters can be deduced, which may lead, in further consequence, to a better understanding and predictability of the growth process.

Erythrocytes, also known as red blood cells (RBCs), are the most abundant type of cells in the human body. Their distribution in the microcirculation determines the oxygen delivery and solute transport to tissues. This process relies on the partitioning of RBCs at successive bifurcations throughout the microvascular network and it is known since the last century that RBCs partition disproportionately to the fractional blood flow rate, therefore leading to heterogeneity of the hematocrit (i.e. volume fraction of RBCs in blood) in micro-vessels. Usually, downstream of a microvascular bifurcation, the vessel branch with a higher fraction of blood flow receives an even higher fraction of RBC flux. However, both temporal and time-average deviations from this phase-separation law have been observed in recent works.

Amongst others, the specific shape and high flexibility of RBCs might cause them to linger, i.e. RBCs can temporarily reside near the bifurcation apex with diminished velocity. This is of critical importance for the cells which are travelling in the capillaries of our circulatory system, whose cross section is sometimes smaller than the erythrocytes diameter. In this presentation, we quantify how the microscopic behavior of RBC influences their partitioning, through combined in vivo experiments and in silico simulations. In vivo experiments are performed in dorsal skinfold chambers of Golden Syrian hamsters, while in silico simulations reproduced the in vivo geometries and are performed using an immersed-boundary-lattice-Boltzmann method.

We developed an approach to quantify the cell lingering at highly-confined capillary-level bifurcations and demonstrate that it correlates with deviations of the phase-separation process from established empirical predictions by Pries et al. More accurately, we demonstrate a linear correlation between the lingering amplitude and deviations from the classical model of erythrocytes partitioning from the literature, both in vivo and in silico. We also show that the lingering is more pronounced in bifurcations presenting a higher curvature at their stagnation point, i.e. the intersection of the bifurcation’s wall and the flow divider plane. This explains why some geometries present stronger lingering than others. Furthermore, we shed light on how the cell membrane rigidity can affect the lingering behavior of RBCs, e.g. rigid cells tend to linger less than softer ones. Taken together, RBC lingering is an important mechanism that should be considered when studying how abnormal RBC rigidity in diseases such as malaria and sickle-cell disease could hinder the microcirculatory blood flow or how the vascular networks are altered under pathological conditions (e.g. thrombosis, tumors, aneurysm).

Blood is a vital fluid that plays an important role in transporting necessary substances and metabolic products for the cells. It consists of about 55% plasma and 45% blood cells where red blood cell (RBC) is the most abundant. Aggregation of RBCs is reversible under normal conditions, but abnormal clusters of RBCs appear in some pathological situations, such as diabetes. The formation mechanism of these clusters is controversial nowadays, but it has been demonstrated that some high-molecular-weight molecules such as fibrinogen and dextran could promote this aggregation.

RBCs surface is covered by glycolipids and glycoproteins, which are collectively called glycocalyx. The alteration of glycocalyx can affect the function of RBCs and it has been reported that a correlation between the enhanced aggregation and the glycoprotein layer (glycocalyx) degradation of the RBCs membrane could be found in diabetes. Other research has also shown that the glycocalyx alteration due to enzyme activities enhanced RBCs aggregation under static conditions by using rat blood. Against this background, glycocalyx alternation is an important point to be studied in this research, especially how glycocalyx alteration can affect RBCs aggregation under flow.

Amylase is a digestive enzyme that plays a crucial role in breaking down polysaccharides. As RBCs are covered by glycocalyx, amylase definitely interacts with RBCs’ glycocalyx. We used pancreatic amylase at various concentrations (0 – 2000U/L) to cleave glycocalyx, mimicking the pathologies of RBCs. We first confirmed the evolution of glycocalyx brush density as a function of amylase concentration with the help of Alexa488-conjugated wheat germ agglutinin (a specific marker of glycocalyx) and confocal microscopy. The results in absence of flow showed that at high concentrations of amylase, the RBCs glycocalyx underwent degradation. Then the aggregation size and morphology were then studied with an inverted microscope equipped with a high-speed camera and a blue light source. These static situation results proved not only the number of RBCs per aggregate was affected by the degradation of glycocalyx but also their morphology was influenced at the concentration of 150kDa dextran in Phosphate Buffered Saline (PBS). The aggregation activities of RBCs under flow were later studied in an artificial polydimethylsiloxane (PDMS) microfluidic circuit with variable pressure difference. The results also showed that the amylase-treated RBCs had significantly more aggregated than non-treated RBCs in 150kDa dextran (15mg/mL in PBS). Cleavage of glycocalyx considerably affected RBCs clusters' size, morphologies, and stabilities under flow.

The mechanism by which amylase caused glycocalyx alteration induces RBCs aggregation is not yet fully understood, and further research is needed to elucidate the exact pathways involved. However, our findings provide valuable insights into the role of glycocalyx alteration in RBC aggregation in absence of flow or in presence of flow and also suggest that glycocalyx degradation by pancreatic amylase has a clear impact on aggregation.

Introduction

Red blood cells (RBCs) are the most abundant cell type in mammals. They are characterized by a high surface-to-volume ratio which gives them very high deformability. When immersed into the flow of blood plasma, this high deformability creates a fascinatingly rich and intricate dynamics due to the two-way coupling between cellular and fluid mechanics. Besides its physiological relevance, e.g. for oxygen delivery, observations of red blood cell shapes have recently been suggested as a tool for medical diagnosis [1, 2]. With some exceptions, typical experimental setups are however limited to 2D microscopic imaging which can represent a serious limitation given the complexity of 3D RBC shapes in blood flow. If properly validated, computational simulations can resolve this issue and provide a deeper understanding of RBC dynamics. Most simulation models which are currently in use, however, are based on experimental validations using static (e.g. optical tweezers) or steady-state properties. In this talk, I will present computational modeling of single RBCs in microchannels focussing specifically on dynamics and transient situations.

Methods

We use boundary-integral computer simulations for the flow computations in the Stokes regime. Membrane dynamics are computed using Skalak and Helfrich laws for shear/area elasticity and bending, respectively. Our models furthermore include the viscous dissipation inside the membrane.

Results

The two prototypical states of a red blood cell flowing in a rectangular microchannel are the croissant and the slipper state [3]. The former is a stationary shape flowing in the center of the channel. The latter is characterized by an off-center flow position and, most importantly, by a constant rotation of the membrane around the liquid interior.

In the first part of the talk, I will briefly introduce computer simulations together with experimental data from our collaboration partners demonstrating that the croissant/slipper dynamics is in fact bistable and that the final RBC shape depends on the initial position. The main focus will then be to develop computational simulations that can quantitatively reproduce the rotation frequency of the membrane in the slipper state. In the second part, I will focus on transient situations such as transitions from croissants to slippers [4] as well as RBCs passing through channels with varying geometries.

I will show that, in order to obtain a computational model which is consistent with this highly diverse range of experiments, it is essential to (i) include membrane viscosity and (ii) employ a higher-than-usual ratio of the RBC interior to that of the outer fluid. Our simulations thus allow a direct and quantitative estimation of RBC membrane viscosity.

[1] Kubánková, M. et al. Biophys J 120, 2838–2847 (2021)

[2] Recktenwald, S. M. et al. eLife 11, e81316 (2022)

[3] Guckenberger, A. et al. Soft Matter 14, 2032–2043 (2018)

[4] Recktenwald, S. M. et al. Biophys J 121, 23–36 (2021)

Sickle Cell Disease (SCD) is a highly prevalent and handicapping genetic disease characterized by hemolytic anemia and unpredictable vaso-occlusive crises (VOC). It is caused by sickle hemoglobin, which, upon deoxygenation and dehydration, polymerizes into fibers in the cytoplasm of red blood cells (RBCs), increasing their cytoplasmic viscosity and remodeling their cytoskeleton. These changes combine to decrease the RBC deformability, a factor involved in the pathophysiology of SCD, with a key role in clinical outcome and occurrence of VOC. Even with reoxygenation, sickle hemoglobin polymerization is only partially reversible. Therefore, the RBC deformability of a sickle patient is highly dispersed as circulating RBCs have undergone a variable number of desoxygenation /reoxygenation cycles. To date, the assessment of RBC deformability require sophisticated devices. A simple method, requiring only a few blood microliters, easy to implement, fast, sensitive to the different parameters governing RBC deformability, and valid for heterogeneous RBC populations is lacking.

We exploit the fact that the individual movement of an RBC under shear flow is an indicator of its deformability. By increasing the flow shear rate, the motion of an RBC suspended in a medium of 40 mPa.s viscosity gradually changes from a rigid-like flip-flopping to a wheel-like rolling motion before transiting, above a critical shear rate set by the cell deformability, to a stationary droplet-like tank-treading motion. The more deformable the cell, the lower the critical shear rate. A given shear rate defines an RBC deformability threshold above which the cell is deformable enough to tanktread. We measured the fraction, f_TT, of tanktreading RBCs in a blood sample, which is the fraction of cells having a deformability higher than this threshold. The lower the applied shear rate, the more severe the deformability test and the fewer RBCs in tanktreading motion is observed in a sample; f_TT is a marker of deformability of the whole RBC sample, which can be tuned by the applied shear rate. We show that this marker is lower in SCD than in control samples. Furthermore, we show that in vitro, this parameter is sensitive to RBC density and hydration, two common features of SCD. We performed a clinical trial on 21 SCD patients by performing weekly repeated f_TT measures over a period of 7 to 23 weeks. We show that f_TT correlates with quantitative biological markers such as level of fetal hemoglobin, reticulocytes count, and plasma LDH level in SCD patients, thus suggesting that a higher hemolysis rate is associated with a lower RBC deformability; f_TT significantly decreases when the patient undergoes a high level of pain, and, surprisingly is increased when patients have an antihypertensor treatment. Moreover, f_TT significantly decreases before VOC onset and arrival at the hospital, then gradually increases during resolution of the VOC. Overall, f_TT could be useful for clinical monitoring of SCD patient to assess treatment efficacy or predict VOC. f_TT is also a good candidate to be an effective marker for diseases associated with RBC stiffening such as malaria.

Over the past 20 years high resolution micro-CT has revolutionized the analysis of trabecular and cortical bone microstructure, facilitating efficient three-dimensional (3D) analysis. Longitudinal, in vivo, applications of this technology, which open up a fourth dimension (4D), in preclinical animal models have primarily focused on trabecular bone, with relatively few studies targeting cortical microstructural features such as porosity. This is despite the growing recognition of the role of cortical bone loss in bone fragility. Two primary factors have contributed to the dearth of in vivo studies of cortical bone: 1) radiation dose, which scales non-linearly with resolution for X-ray computed tomography, is considerably higher when targeting cortical pores as they are smaller than trabeculae; and 2) small rodents, including mice and rats, exhibit little, if any, intra-cortical remodeling comparable to larger vertebrates, including humans. This presentation details progress over the past decade to overcome these hurdles and advance in vivo imaging of cortical bone porosity. Specifically, we have pursued 4D tracking of the remodeling spaces that lie at the heart of basic multicellular units (BMUs) to shed new light on the spatio-temporal regulation of remodeling. To overcome the challenge arising from radiation dose, our group has developed synchrotron-based imaging protocols at the Biomedical Imaging and Therapy facility of the Canadian Light Source synchrotron which capitalize on in-line phase contrast micro-CT. This approach has enabled enhanced detection of cortical bone porosity, while minimizing radiation dose (1-5 Gy) and scan time (<1 minute). To address the limitation of small rodents, we have developed rabbit models of elevated cortical bone remodeling/porosity, including ovariohysterectomy, parathyroid hormone (PTH) dosing and glucocorticoid dosing. By combining synchrotron-based imaging with rabbit models, we have achieved the first 4D tracking of remodeling spaces and direct analysis of their rate of advance (linear erosion rate; LER). We have discovered that LER is reduced in animals undergoing PTH dosing compared to those withdrawn from PTH, a potential mechanism by which coupling of resorption and formation are enhanced by this hormone. We are actively developing novel 3D morphological tools for assessing remodeling spaces, to measure the extent and timing of the resorptive, reversal and formative phases of BMUs and to examine potential impacts of PTH and other interventions. We are also exploring these morphological approaches as a means to facilitate comparative analysis between rabbits and humans to establish the validity (and limitations) of this experimental platform for the study of osteoporosis etiology and treatment. Finally, a key objective of our team is the integration of 4D data into computational models of cortical bone remodeling to accelerate discovery related to the improvement of osteoporosis treatment. Direct study of the spatio-temporal regulation of BMUs is enabling a shift towards the assessment of remodeling at its fundamental level, rather than making inferences from accumulated microstructural changes.

Bone is a defining feature of all vertebrates, present in more than ~65,000 extant species today. Even cartilaginous fishes, such as sharks and rays, and thought to have possessed a bony skeleton early on in their evolutionary history before secondarily losing it. For palaeontologists studying the history of vertebrates, bone is the most commonly preserved trace of an animal and provides the literal and figurative backbone upon which to reconstruct past life.

The evolution of the earliest four-limbed terrestrial animals from water-dwelling fishes is hailed as one of the greatest ‘steps’ in evolution and is thought to have occurred during the Devonian Period, more than 360 million years ago. Significant adaptations were required in the bodies in the first tetrapods to enable them to leave the water and conquer land, including the appearance of limbs, digits and lungs, and concomitant strengthening of the limb bones, girdles and axial skeleton.

Recently, several novel imaging and computational methods such as synchrotron and computed tomography (CT) have transformed palaeontology enabling non-destructive examination of rare and irreplaceable fossils down to histological scale, the creation of digital models for complex biomechanical analyses, and revealing hitherto unseen internal anatomy. In this talk I will detail several recent research projects examining adaptations and remodelling in the bony skeletons of fish and early tetrapods as revealed by novel imaging methods.

First, high-energy computed tomography reveals the skeleton of the pectoral fin in Elpistostege, a tetrapod-like fish from the Late Devonian of Canada. The pectoral fin endoskeleton contains a humerus, radius, ulna, ulnare and four proximodistal rows of radials and two distal rows organized as digits. This tetrapod-like pattern is retained within a fin with distal lepidotrichia but represents an important stage in the early evolution of the vertebrate hand.

Secondly, the humeral microarchitecture of stem-tetrapods, batrachians, and amniotes was examined using three-dimensional synchrotron virtual histology. We show that a centralised marrow organisation (to enable haematopoiesis, the production of blood cells, as exhibited in living amniotes) did not arise during the water-land transition as originally hypothesised, but arose considerably later in Permian amniotes.

Lastly, Finite Element Analysis (FEA), a technique common in engineering and mathematical modelling to predict how structures react to various forces, is herein applied to address questions about form-function relationships in the evolution of vertebrates. FEA is applied to several of the lungfishes from the Late Devonian Gogo Formation, Australia, the most diverse lungfish assemblage in the world. The >10 described species display extreme variation in skull, dentition and mandible morphology which has been proposed as a driver for their success. Here using FEA we show that robust forms exhibit higher strain tolerances during feeding compared to more gracile forms. Our results demonstrate that biomechanical function and feeding performance are constrained by mandible morphology in Devonian lungfish.

These three case studies show significant bone adaptations observed during early fish-tetrapod evolution (across histological to gross morphological scale) as revealed by new tomographic computational methods, and highlight the universality of bone as a vertebrate tissue.

John Clement was my student at UCL. When he graduated and wanted to undertake further training in pathology, I advised him of a vacant position at The London Hospital Medical College [later QMUL], which he took, and thus entered micro-anatomical and forensic research. When he removed to Melbourne, he established the renowned femur collection which enabled many international collaborators to access first class, well documented human bone research material. At UCL, we established a facility for determining the fabric level mineralisation density of bone at the sub-micron scale using quantitative backscattered electron imaging in an automated digital scanning electron microscopy system [qBSE-SEM]. We also introduced confocal scanning optical microscopy [CSLM] to bone studies and developed methods for marrying the information contents from these important methodologies. At the same time, Jim Elliott developed X-ray microtomography [XMT] at QMUL and we correlated all these methods.

John provided us several femur samples, which like most other research material at the time we embedded in PMMA, en route to smart imaging with SEM and CSLM. John, with David Thomas, took to XMT in Melbourne. Unable to compete with XMT for large volume 3D imaging for the vascular space compartment in bone, we took to casting aside the bone in our precious blocks by dissolving it with sequential treatments with HCl and NaOCl solutions, leaving a cast of the marrow spaces and the osteocyte lacunar-canalicular space. The latter was so abundant that it blocked visibility of the blood vessel canal spaces, so we destroyed it by ultrasonication to clean the residual cast. This was in earlier days coated with gold to give surface conductivity and an enhanced BSE signal, but latterly we have simply used uncoated samples at say 50Pa chamber pressure to circumvent charging problems.

The resin casts allow us to image the space compartment in bone at a resolution far superior to that obtainable by XMT, allowing new insights into growth, modelling and remodelling processes in compact bone tissue.

The recent decade our understand of the cellular organization of human bone remodeling have been revised, providing a new perspective on this process during physiological and pathophysiological conditions, and the critical remodeling steps contributing to bone loss. We know now that bone remodeling comprises three successive phases: (1) a short initial resorption phase by primary osteoclasts, (2) a longer reversal-resorption phase with intermixed reversal cells (osteoprogenitors) and secondary osteoclasts, and (3) a subsequent formation phase with mature bone-forming osteoblasts. Moreover, it is well established that trabecular bone remodeling events are separated from marrow cavity by canopy cells (osteoprogenitors), which is part of the bone marrow envelope enwrapping the marrow cavity.

In this study, we focus on the close interplay between osteoclasts and osteoprogenitors (bone lining cells, reversal cells and canopy cells), which play a critical role in coupling of bone resorption to the subsequent bone formation. During the reversal-resorption phase, the osteoprogenitors colonizing eroded surfaces gradually differentiate into bone-forming osteoblasts and expand to a critical density required for initiation of bone formation. Osteoclasts are in direct cell-cell contact with these osteoprogenitors (involving Semaphorin 4D - plexin B1 interaction) and secrete coupling factors that promote these osteoprogenitors (involving clastokines like LIF, PDGFB, TRAcP5b), which on the other hand express osteoclast-promoting RANKL and degrade the collagen debris left behind by the osteoclasts. These osteoprogenitors also interact with the mechanosensing osteocytes, and there is evidence supporting that mechanical conditions likely sensed by osteocytes may affect the interplay between osteoclasts and osteoprogenitors. Indeed, osteoclasts, osteoprogenitors and osteocytes form an overlook partnership, which is critical to the coupling of bone resorption and formation within each bone-remodeling event.

Tissue specific metabolic rates for various organs and organ systems are typically calculated for multiple organ tissues as given constants of kcal/kg expended per day. These constants (Kᵢ) are multiplied by an organ’s mass (Tᵢ) and add up to the whole body resting energy expenditure: REE = Σ(Kᵢ × Tᵢ).

Organ and tissue mass-specific metabolic rates are measured using a variety of techniques such as magnetic resonance imaging, calorimetry, or by a respirometer. Within the category of ‘residual mass’ in metabolic rate calculations often includes the skeleton, blood, skin, and other small tissues and organs. Problematically, skeletal tissue is lumped in with other organs and tissues for the metabolic rate constant of residual mass since it is not considered an expensive tissue with a high metabolic rate. But this ignores the exponential scaling of osteocyte density with body mass reported in mammals.

Metabolic rate scales at the ¾ power of body mass. This same scaling can be interpreted as a ratio of the average metabolic rate of a cell, Bc, to the average cell size, mc: B/M = Bc/mc. The average cellular metabolic rate to average cell size decreases as body size increases so that as body mass changes, there is a tradeoff between cell size and cell metabolic rate. However, osteocytes generally fall within the narrow range of 5-20 µm diameters, but comparative mammalian osteocyte density increases nonlinearly with body mass. Thus the tissue specific metabolic rate of bone must vary enormously across the Class.

Bone, unlike other organ tissues, is unique in that roughly 95% of its cells are trapped in an inorganic mineral matrix at densities per unit volume that depend upon growth rate, body mass, and thus also life history. In our intraspecific research on human bone, osteocyte density scales positively with body height, indicating that larger individuals have higher osteocyte densities for orchestrating bone formation and mineralization, which thus also have higher metabolic rates. It appears that the same increase in energetic efficiency observed in interspecific comparisons of the mass specific metabolic rate of bone at larger body sizes also characterizes body size categories among humans. Long period biological rhythms that regulate rates of cell proliferation explain some aspects of human body size variability.

Retrieving metabolic rates or rates of oxygen consumption from living organisms is a difficult task that has not been fully explored in bone. Our presentation will examine the scaling of osteocyte density and its implications for the mass-specific metabolic rate of bone.

Bone is a hierarchical, composite connective tissue that can adapt and change through the lifespan of an individual. Its mineralized matrix holds within it information about life history, functional adaptation, health and disease, leaving tell-tale signatures (e.g. in terms of tissue composition, distribution and orientation and the spaces within it) as a result of the dynamic processes of modeling and remodeling. Contributions from fields as varied as paleontology, forensic science, anthropology, orthopedics, materials science and engineering, combined with advances in imaging and analysis methodologies have resulted in a deeper understanding of this complex internal structure of bone and its relationship to attaining and maintaining bone quality and responsiveness to mechanical loading through the lifespan.

This presentation reviews multidisciplinary research that grew from collaborations with Dr. John Clement and the Melbourne Femur Collection, and expanded into other anatomical locations and samples with the aim of (1) quantifying patterned organization in cortical bone microstructure and properties through ontogeny and with aging using correlative light and electron microscopy as well as 3D techniques such as microCT, (2) relating this organization to the underlying processes of bone modeling, remodeling and establishment of bone morphology, and (3) applying these techniques to bioarchaeological and clinical contexts. The presentation will highlight how new methods of visualization at different length scales can enhance not only this kind of basic science research, but also be used in the education of future researchers and clinicians.

The fetal membranes are an important biological structure for pregnancy that surrounds and protects the fetus. They are layered structures, comprising the amnion, the chorion, and part of the maternal layer decidua. During gestation, they undergo complex microstructural changes, such as the weakening of the tissue in preparation for delivery. Little is known about the influence of the collagen fibers organization within the amnion in the rupture process of the fetal membranes. Non-crystalline X-ray diffraction (XRD) has been applied to study the collagen organization in some soft tissues, such as the cornea, breast, bone, cartilage, and, more recently, amnion. Some studies have suggested that the loss of the regular structure of the collagen fibers in the amnion may represent a contributing factor to PPROM. This work aims to analyze the biomechanics of fetal membranes under different off-plane collagen fibers within it through the potential of numerical analysis to understand whether it has a potential impact on preterm pre-labor rupture of the fetal membrane (PPROM).

A multilayer model of the fetal membrane was developed based on a robust inflation mechanical test dataset. In terms of constitutive models, the amnion was characterized by the modified version of the Buerzle-Mazza constitutive model (μo=2.4MPa, q=2.96, m5=0.463, m2=0.00228, m3=41.12, m4=1.27, N=32). The chorion (E=1MPa, υ=0.41) and the decidua (E=1MPa, υ=0.49) were characterized by elastic linear properties.

The evolution of the maximum principal stress curves with pressure in the amnion when the off-plane angle is set to 0° or 10° is different from the cases where the same parameter is set to 20° or 30°: for the first two values, the curves always increase throughout the entire simulation; for the 20° and 30° angles, the stress is 0 for smaller pressures and a rapid increase in stress is verified for higher pressures. The maximum principal stress is larger when the angle is changed from 0° to 10°, and from 30° to 20°. In the chorion layer, the maximum principal stress at the apex of the membrane increases with the off-plane angle.

Different off-plane fiber angles had a strong impact in terms of maximum principal stress in both layers, especially in the mechanical dominant layer amnion. In terms of PPROM, it is very likely that certain off-plane collagen fibers that lead to higher maximum principal stresses potentiate the rupture of the membrane. These results highlight the potential of our model to characterize the biomechanics of fetal membranes under different physiological conditions.

Pelvic floor disorders, including uterine prolapse and urinary incontinence, have a significant impact on women’s quality of life and well-being, affecting between one-third and one-half of all women [1]. For uterine prolapse, risk factors include advancing age, obesity, physical inactivity, and parity. The disease is caused by a weakening of pelvic floor musculature and other supporting anatomical structures, such as the uterine ligaments. Currently, tissue changes leading to the development and progression of uterine prolapse through mechanical and molecular alterations are understudied. Since aging and hormonal changes may affect components of the extracellular matrix, including collagen fibrils and their nanostructure (d-spacing) [2], it is crucial to investigate changes in the collagen backbone of the uterine support structures to gain a deeper understanding of the development and progression of uterine prolapse. As part of the uterine support structures, the broad ligament connects the uterus to the lateral pelvic walls and is often described as a double layer of peritoneum. Consequently, our aim was to investigate the characteristics and nanoscale behavior under loading of collagen in the broad ligament.

Broad ligament samples were obtained according to local ethics regulations from women of different ages (66 ± 22 years) post mortem. Quadratic samples (15 × 15 mm²) were subjected to biaxial stretching with simultaneous microfocussed ultra-small-angle X-ray scattering (USAXS) to determine changes in collagen fiber orientation and d-spacing with deformation (MiNaXS beamline P03/PETRA III, DESY) [3]. Deformation was increased in a stepwise manner (0%, 5%, 10%, 15%, 20%). For each step, a map of 1 × 1mm² with 25 individual USAXS measurements was created. In addition, samples adjacent to the USAXS-samples were prepared for histological assessment of collagen orientation, elastin and proteoglycan content, and cellular properties. Slices were cut in frontal, sagittal, and horizontal orientation to allow collagen fiber assessment from different directions.

Radial integration of scattering data indicated the presence of two main orthogonal collagen fibril orientations. With increasing tissue strain, collagen d-spacing (an indicator of fibril strain) and fibril alignment increased while fibril thickness decreased (all p<0.05). Preliminary histological evaluation from slides with orthogonally oriented cuts confirmed the presence of multiple main collagen fiber orientations. Collagen fibers exhibited crimping. Further, the broad ligament samples were cell rich and included small vessels.

Histological evaluation confirmed the presence of multiple predominant fiber orientations as indicated by USAXS. This is in agreement with descriptions of fiber distribution in peritoneum [4]. Compared to reported fibril strains in tendons [5], the observed fibril strain in the broad ligament was smaller. This could be explained by the much higher alignment of collagen fibers to a unidirectional loading axis in tendons. The presented findings provide a more detailed understanding of collagen characteristics of the broad ligament and may contribute to the development of biomechanical models of the uterine support system.

References:

[1] MacLennan et al., 2000, DOI: 10.1111/j.1471-0528.2000.tb11669.x.

[2] Fang et al., 2012, DOI: 10.1038/jid.2012.47.

[3] Euchler et al., 2022, DOI: 10.1088/1742-6596/2380/1/012109.

[4] Liu et al., 2017, DOI: 10.1016/j.biomaterials.2016.11.041.

[5] Barreto et al., 2023, DOI: 10.1016/j.matbio.2022.11.006.

Childbirth trauma during the second stage of labor is a prevalent concern, affecting millions of women worldwide. Levator ani muscle (LAM) trauma, which can impact 6-40% of women undergoing vaginal delivery, particularly nulliparous individuals, is a prevalent injury arising from childbirth. LAM trauma can lead to persistent morbidity and the necessity of future surgical intervention for 10-20% of patients. However, predicting and diagnosing pelvic floor lesions can be challenging. To address this issue, biomechanical simulations have been used as a helpful tool to evaluate pelvic floor muscle (PFM) injuries. However, using the finite element method (FEM) for these simulations can be a time-consuming process. Thus, it is important to explore alternative techniques for applying these methods in a clinical setting. One promising approach is to use machine learning (ML) algorithms, which can leverage simulation data to offer faster results. The present study aims to develop a ML framework to predict stresses on the PFM during childbirth by training ML algorithms on FEM simulation data.

To generate the data, 2744 childbirth simulations were performed in Abaqus software, in which the material parameters of the constitutive model used to characterize the PFM were changed. The constitutive parameters varied between ranges of values ​​according to the literature, thus allowing to characterize the pelvic floor of most women. A total of 1715 simulations were successfully completed.

A dataset was created in which each observation corresponds to a node of the pelvic floor in one simulation. Specifically, 46 nodes located in the inferior portion of the PFM were selected, which is the region that undergoes the most stretching during childbirth. Features such as node number and position, initial coordinates, and material parameters were used for training. Four ML models, namely Random Forest (RF), Extreme Gradient Boosting (XGBT), Support Vector Regression (SVR), and Artificial Neural Networks (ANN), were chosen for the study. A training and test set were created with a 90/10 split, recurring to a stratification method to guarantee the same feature distribution in both sets. Subsequently, hyperparameter optimization with cross-validation was performed. The models' performance was measured by the mean squared error (MSE) and the mean absolute error (MAE).

The stress values were measured at the moment of maximum stretch, and ranged from 0 to 30 MPa. The results for predicting the maximum principal stresses showed that the ANN produced the best outcomes, with a MSE of 0.112 and a MAE of 0.191. Conversely, the SVR model had the highest error, with a MSE of 0.444 and a MAE of 0.356. Both tree-based algorithms performed reasonably well and were closer to the outcomes achieved by the ANN. The ANN is capable of making predictions in approximately 120 milliseconds, indicating its potential for real-time applications.

The current work represents an advance in the field of childbirth computational simulations using artificial intelligence tools. The ability to predict the stresses suffered by the woman on the pelvic floor immediately before or during childbirth could aid in medical decision-making and in the identification of non-visible injuries.

Tissue development in normal and pathological conditions is driven by a set of biological phenomena. Growth is characterized by a change of mass which might be positive (tissue growth) by cell division, cell enlargement, and extracellular matrix secretion, or negative (tissue atrophy) by cell death, cell shrinkage or resorption [1]. Cancer is a malignant pathology characterized by accelerated and uncontrolled cellular growth and proliferation [2]. In 2020, it was estimated 19.3 million newly diagnosed cases and almost 10.0 million cancer deaths, worldwide [3]. The development of solid tumors is a multi-factor biological process, as it is influenced by molecular and genetic factors, cell-cell and cell-extracellular interactions, and vascularization (in particular oxygen and nutrients supply) [4].

Several experimental and theoretical efforts have been taken to unravel the mechanisms of tumor growth. Constitutive mechanics models are one of these approaches as they can help describe mass changes and stress development during the events of tumor progression [4]. From a biomechanical perspective, solid tumors are hyperelastic, compressible, anisotropic materials with a mechanical behavior both space and time-dependent, whose mass varies over time [5]. In this work, the goal is to establish a computational framework to model tissue growth to be latter applied for solid tumors.

Initially, a constitutive model to describe the anisotropic growth of a solid mass is implemented [6] considering the multiplicative decomposition of the deformation gradient into an elastic and a growth contribution [7]. In this first attempt, anisotropic growth is considered a strain-driven process with a privileged direction for growth due to the presence of a micro-structure [8]. To ease the use of this model in more complex scenarios, the constitutive model is implemented into the finite element software Abaqus® using a user-defined subroutine (UMAT). To validate the UMAT subroutine, the numerical solution for a single unitary hexahedral finite element is computed for a set of deformation cases in the software ABAQUS® and compared to the analytical solution [6].

Finally, the UMAT is implemented in an in-silico model of a solid tumor surrounded by the extracellular matrix (simplified by a parallelepiped) and cyclic uniaxial and biaxial stretch scenarios are applied. The Cauchy stresses are recorded in the direction of the applied displacements as function of the stretch. At the unload final state, the embedded tumor presents a positive volume change as the tumor grows irreversibly in an anisotropic manner. Future steps include the incorporation of a tumor growth law derived by experimental evidence and the implementation of a stress-driven anisotropic growth evolution.

References

[1] Ambrosi D, Guana F. Stress-Modulated Growth. Mathematics and Mechanics of Solids 2007;12:319–42. https://doi.org/10.1177/1081286505059739.

[2] Ambrosi D, Amar M Ben, Cyron CJ, DeSimone A, Goriely A, Humphrey JD, et al. Growth and remodelling of living tissues: Perspectives, challenges and opportunities. J R Soc Interface 2019;16. https://doi.org/10.1098/rsif.2019.0233.

[3] Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209–49. https://doi.org/https://doi.org/10.3322/caac.21660.

[4] Xue S-L, Li B, Feng X-Q, Gao H. Biochemomechanical poroelastic theory of avascular tumor growth. J Mech Phys Solids 2016;94:409–32. https://doi.org/10.1016/j.jmps.2016.05.011.

[5] Ramírez-Torres A, Rodríguez-Ramos R, Merodio J, Penta R, Bravo-Castillero J, Guinovart-Díaz R, et al. The influence of anisotropic growth and geometry on the stress of solid tumors. Int J Eng Sci 2017;119:40–9. https://doi.org/10.1016/j.ijengsci.2017.06.011.

[6] Vila Pouca MCP, Areias P, Göktepe S, Ashton-Miller JA, Natal Jorge RM, Parente MPL. Modeling permanent deformation during low-cycle fatigue: Application to the pelvic floor muscles during labor. J Mech Phys Solids 2022;164:104908. https://doi.org/10.1016/j.jmps.2022.104908.

[7] Rodriguez EK, Hoger A, McCulloch AD. Stress-dependent finite growth in soft elastic tissues. J Biomech 1994;27:455–67. https://doi.org/10.1016/0021-9290(94)90021-3.

[8] Soleimani M, Muthyala N, Marino M, Wriggers P. A novel stress-induced anisotropic growth model driven by nutrient diffusion: Theory, FEM implementation and applications in bio-mechanical problems. J Mech Phys Solids 2020;144:104097. https://doi.org/10.1016/j.jmps.2020.104097.

Acknowledgements: The authors gratefully acknowledge the support from the Portuguese Foundation of Science under the grant SFRH/BD/09480/2022 and the funding of the research project PTDC/EME-APL/1342/2020.

The Finite Element Method (FEM) is a powerful tool that enables the simulation of many complex engineering problems. In complex analysis, such as in the modelling of biological and biomechanical phenomena, it is necessary to specify the constitutive equations that describe the biomechanical behaviour of such materials.

Soft tissues and other biological materials subjected to large deformations present an extremely nonlinear behaviour, which makes the constitutive modelling of such materials a complex task and expensive in terms of time and computational resources. Alternatively, it is possible to use surrogate models, which consist of models that replace the traditional and expensive main models to overcome some computational limitations.

Such surrogates can learn the behaviour of the soft tissue when trained on previously acquired data and then replace the expensive numerical models. To develop the surrogates, Artificial Neural Networks (ANNs) will be trained from a large dataset with the deformations (inputs) and the corresponding stresses (outputs). Then, the weights and biases of the trained model will be used to write the forward pass equations in Fortran to implement a general user material subroutine (UMAT) for the Finite Element software ABAQUS. The surrogate models will be tested under homogeneous deformation cases and under more complex examples, where it is shown the values of the maximum principal stress obtained with a conventional UMAT and with the ANN. In the future, this data driven approach is going to be applied to soft tissues, such as the pelvic floor muscles.

The female pelvic cavity is a very complex anatomical region. Pelvic dysfunction, especially pelvic organ prolapse (POP) and urinary incontinence (UI), have a negative impact on women's lives, and it happens when the support mechanisms of the pelvic cavity become fragile. UI has a prevalence of a up to 28%, with stress urinary incontinence (SUI) being the most common form [1,2], characterized by involuntary urinary leakage during physical strain, coughing or an increase in intra-abdominal pressure (IAP). Existing treatments for these disorders are divided into conservative and invasive. The last ones consist of surgical interventions and should be used in patients in whom the first treatments did not work, or when the severity of the dysfunction is high.

SUI occurs when the intravesical pressure exceeds urethral resistance at which the urethra has the ability to remain closed [2]. Furthermore, it comes to a point when neither the pubourethral ligaments (PULs) [3] and the arcus tendineous fasciae pelvis (ATFP) [4] can stabilize the bladder neck (BN) [5]. Assessment of BN mobility in patients with SUI is essentially clinical, however, the imaging techniques such as ultrasound (US) and magnetic resonance imaging (MRI) are used as a method for evaluating this characteristic. The outcomes of radiographic images have been crucial and used as input for numerical methods.

The aim of the present study was to establish the IAP values and the in vivo biomechanical properties of the bladder tissue for two distinct groups (continent women and women with SUI). The numerical simulations of Valsalva maneuver were performed, applying the Ogden hyperelastic constitutive model to the bladder and also the inverse finite element analysis (FEA).

This study focuses on adapting an inverse FEA to estimate the in vivo properties of the bladder, using a constitutive model of the female pelvic cavity and MR images acquired at rest and during the Valsalva maneuver, for two distinct groups (continent and incontinent women). The bladder neck’s displacements were compared between computational simulation and MR images.

The results of the FEA showed that the bladder tissue of incontinent women have the highest stiffness approximately 47% higher when compared to continent women.

References

1. Jansson MH, et al. Stress and urgency urinary incontinence one year after a first birth—prevalence and risk factors. A prospective cohort study. Acta Obstetricia et Gynecologica Scandinavica 2021; 100(12):2193–2201.

2. Falah-Hassani K, et al. The pathophysiology of stress urinary incontinence: a systematic review and meta-analysis. International Urogynecology Journal 2021; 32(3):501–552.

3. Kefer JC, et al. Pubo-urethral ligament transection causes stress urinary incontinence in the female rat: A novel animal model of stress urinary incontinence. Journal of Urology 2008; 179(2):775–778.

4. Iyer J, et al. Introduction and Epidemiology of Pelvic Floor Dysfunction. In: Rane A, Rane A, Durrant J, Tamilselvi A, Sandhya G, eds. Ambulatory Urology and Urogynaecology. Wiley Online Library; 2021.

5. Occelli B, et al. Anatomic study of arcus tendineus fasciae pelvis. European Journal of Obstetrics and Gynecology and Reproductive Biology 2001; 97(2):213–219.

Tumour growth is a force-sensitive process, regulated in part by mechanical feedback from surrounding tissue [1]. Such mechano-responsiveness can govern tissue-specific risk and progression of cancer. However, the underlying biomechanisms by which mechanical loading influences cellular growth and proliferation have not yet been uncovered. In this work, we propose a novel computational framework to determine how the feedback between growth, mechanical loading, and cell-cell exchange flow could restrict tumour cell proliferation.

Non-cycling cell volume is regulated by an interplay between energy-consuming pumps, mechanosensitive ion channels, and actomyosin tension that coordinate to manipulate cellular osmotic pressure. Growth may be induced by impermeable solute synthesis and associated increases in osmotic pressure, which drives water intake. Through a combination of novel finite element analyses (FEA) and agent-based modelling (ABM), we investigate how internal and external stresses arising from active cell tension and mechanical loading could oppose limit cell division by restricting growth below a critical mitotic volume.

Our model predictions for osmotic control of division indicate that cell cycle synthesis drives growth, and that compressive loading can limit the potential for a cell to surpass the size checkpoint for division. Simulations further suggest that exchange flow through gap junctions in connected cells can restrict osmotically-regulated cell growth and subsequently restrict proliferation. Using our integrated FEA-ABM modelling framework to characterize multicellular interactions and matrix loading, simulations reveal that increasing matrix stiffness reduces the rate of cell proliferation in tumour spheroids, due to the emergent stress-sensitivity of cell growth and division. Mean cell pressure is predicted to converge to a critical value, independent of matrix stiffness, at which proliferation is inhibited. Cells at the tumour core are revealed to experience higher stress than peripheral cells, as supported by data from excised tumours [2]. Overall, simulations suggest that tumour spheroid size reduces with increasing matrix stiffness in a stress-dependent manner, supported by our experiments of tumour spheroid growth in hydrogel.

Our analyses suggest that stress-dependent tumour growth emerges from a constraint on osmotically regulated cell growth, whereby cells cannot obtain a critical mitotic volume due to external loading. Simulation of multicellular proliferation with coupled finite element and agent-based models provides unique insight into the evolution of macro-scale tissue behavior and mechanosensitive growth, with broad applications to patient-specific cancer diagnosis.

References

1. M. Kalli and T. Stylianopoulos. Front. Oncol. 8 (2018)

2. H.T. Nia, et al. Nat. Biomed. Eng. 1 (2017).

One of the crucial steps of cancer metastasis is the extravasation of cancer cells from the vascular network into the neighbouring organs. The endothelial cells, which form the inner lining of blood vessels, develop homophilic bonds with neighbouring cells with the help of VE-cadherin proteins and form a preventive barrier against extravasation. But there is strong evidence of such cell-cell bonds abruptly broken in cancer patients that helps the trans-endothelial migration of cancer cells [1]. A diseased state leads to the abnormal dynamic behaviour of VE-cadherin junctions resulting in a larger gap size between cells and altering the frequency of their opening, a state termed hyperpermeability. Chemical cues such as Thrombin, along with mechanical properties of the cell and extracellular matrix (ECM) are known to affect permeability by regulating the actin-cytoskeleton machinery of endothelial cells. But details of how the coupled mechano-chemical stimuli affect endothelial permeability are not yet understood completely.

We have developed a continuum-level mathematical model to study the dynamics of such gap formation. The model considers the mechanical equilibrium of the cell by coupling the VE-cadherin to Actin-cytoskeleton. Taking inspiration from contact mechanics, we model VE-cadherin as a cohesive surface with damage following a traction-separation law. Actin cytoskeleton follows the tensegrity principle and is coupled with microtubules following non-linear material laws. Cells are modelled in 3D and hence we can study the behaviour of micro-vessels, which are hard to study experimentally. The ECM-cell interaction also follows a traction-separation law, allowing us to understand the role of ECM stiffness in cell-cell interactions. We use a non-linear finite element solver to solve the mechano-chemical coupled system of equations.

Numerical simulations showed that micro-vessels exhibit higher permeability than their planar monolayer counterparts. In-silico studies show that the permeability of micro-vessels increases with the stiffness of the extracellular matrix. Interestingly, it was observed that shear between cells is responsible for variation in permeability between bi-cellular and tri-cellular junctions, explaining the phenotypic differences observed in experiments. Simulations also show that permeability is higher in those regions of micro-vessels with high shear stress fluctuations, matching the observations that atherosclerotic plaques are usually formed in regions where blood flow is disturbed. The novel mathematical modelling framework is capable of simulating already known results along with developing several testable hypotheses. The versatile model can thus be used in a variety of studies, from the role of VE-cadherin in the extravasation of cancer cells to understanding the effect of nanoparticles in effective drug delivery.

1. T Tomita et al., “Regulation of vascular permeability in cancer metastasis”, Cancer Sci, 2021

Cell mechanics of living cells are vital for many cell functions, including mechanotransduction , migration, and differentiation. It is known that changes in cell mechanics are often correlated with disease progression. Cell-matrix adhesion is important for the patterning, integrity and homeostasis of tissues and may provide a target for therapy, for example in cancer metastasis. Cell mechanics and adhesion between cells and the matrix are also important for tissue engineering. Therefore, it is important to study cell mechanics and cell-to-material adhesion. However, simultaneous characterization of the cell-to-material adhesion and viscoelastic properties of the same cell is challenging. In this study, we present a new approach to simultaneously determine these properties for single cells, using Microfluidics-based Atomic Force Microscopy.

During the approach period, the cantilever touches the cell and rests for a few seconds, where the cell relaxation curve is recorded. Then, negative pressure is applied to the cell through the cantilever, and the hollow cantilever grabs the cell and pulls it away from the materials surfaces. During the relaxation period, the viscoelastic properties of the cell can be determined where the contact area between cell and the cantilever can be visualized and measured. The adhesive force and adhesive energy for cell-to-material can be determined based on the detachment curve when pulling the cell away from the substrate.

To explore the correlation between cell-materials adhesion and cell migration, the single cell tracking method was adopted. ImageJ with TrackMate plug-in was employed to identify the cells and then analyse their path. The mean migration speed is calculated by dividing the displacement of the cell between each frame by the time between frames to obtain the velocity of the cell moving between those two frames, and then averaging the movement velocity over the entire motion trajectory.

To reveal if the mechanical properties of biomaterials may affect the biomechanics of the breast cancer cells (MCF-7 cells), typical biomaterials Polydimethylsiloxane (PDMS) were used. The curing agent-to-base ratio was varied to yield PDMS with different stiffness. These PDMS samples were treated with Ozone to change the surface wetting properties to enable better cell attachment and proliferation.

The elastic modulus for MCF-7 cells on petri dish, determined here, is consistent with what was reported in literature. We have also demonstrated that both elastic cell moduli and viscosity of breast cancer cells (MCF-7 cells) can be affected by the mechanical properties of PDMS with different stiffness. The individual cell, which has stronger adhesion to the materials, appears to be stiffer. The cells seem to have stronger adhesion and slower migration when interfacing with stiffer materials.

According to the World Health Organization (WHO), prostate cancer (PCa) is the second most common cause of cancer worldwide¹.

Accurate risk stratification is crucial in treating prostate cancer. However, distinguishing high-risk from low-risk cases represents a significant challenge. Current diagnostic tools, such as blood-based prostate specific antigen (PSA) tests, digital rectal examination findings, and biopsies, often fail to differentiate the aggressiveness of cancer, leading to misclassification and under/over treatment². Consequently, new tools are urgently needed to predict outcomes and improve clinical decision-making to prevent unnecessary invasive treatments and associated morbidity³.

Therefore, the aim of this study is to develop an in-silico model of PCa tumour growth or degrowth that will predict prostate cancer outcome and provide an improved prognosis.

For the simulation of PCa tumour growth, a multispecies model of partial differential reaction-diffusion equations coupled with the mechanics of continuous media is here utilised. This multispecies model represents the phenomenological behaviour of PCa and its cellular processes, including proliferation, differentiation, and apoptosis. From these processes, the model simulates the evolution of the geometry, the distribution of tumour and healthy cells and stroma in the entire prostate. The model is implemented via the Finite Element Method. Magnetic Resonance Images (MRIs) are available to reconstruct the geometry of the prostate in a patient-specific manner. These images provide in-depth knowledge of the level of cellularity (DWI sequences) and vascularization (DCE sequences) of the tissue, which serve as input data to the model⁴. Real patient data are available for the validation of the model at a later time point.

The model has been tested on two clinical cases, showing promising preliminary results. The distribution of cellularity and PSA evolution are compared. Over time, the cellularity of the prostate obtained from the simulations reproduces closely the clinical values, with a similar spatial distribution. Furthermore, quantification of global cellularity does not give relative errors greater than 7%. On the other hand, the results obtained from the simulation of the PSA evolution suggest a close correlation with the observed clinical evolution, and the model can accurately predict the PSA's increase or decrease. However, further refinement of the model is necessary to improve its predictive capacity.

A patient specific model has been applied to simulate tumour growth in PCa. Moreover, this model can give predictions of the evolution of the PSA, which is the principal biomarker used by clinicians for the PCa follow-up. In future developments, different treatments will be incorporated which will allow for a better prognosis. Finally, the model will be further validated with additional patient-specific clinical data.

Acknowledgements

This publication is part of the project PLEC2021-007709 (ProCanAid), funded by MCIN/AEI/10.13039/501100011033/ and by the European Union NextGenerationEU/PRT and in collaboration with IISLAFE and QUIBIM. SHR was supported by the Government of Aragon (2019-23).

References

[1] Rebello, R.J., et al. Nat Rev Dis Primers. 7, 2021.

[2] Elwenspoek, M. et al. JAMA Netw open. 2:198427-198427, 2019.

[3] Moore, C.M. et al. Eur Urol, 64:544-52, 2013.

[4] Saiz-DeMena, D. et al. Eng Comput 3849–3865, 2022.

Tumour growth is usually hampered by poor blood supply when the tumour reaches a critical size. Under these circumstances, the process of growth tends to stall. When this supply is not verified, tumour cells may begin to secrete angiogenic growth factors that activate the process of angiogenesis, which is characterized by the formation of new blood vessels towards the tumour. Understanding this mechanism is one of the keys to comprehending tumour progression. To help in achieving this objective, computational models to describe angiogenesis in this context have emerged and have shown their potential to aid in this goal. Numerical methods are usually associated with these models, and numerous ones are documented in the literature. One of the most famous is Smoothed Particle Hydrodynamics (SPH).

The objective of this work was to enhance a 3D algorithm created by the authors, which simulates cell proliferation using the SPH algorithm along with the process of angiogenesis. In this algorithm, angiogenesis is triggered by VEGF concentration, and to test it, different concentrations and locations for the focus of concentration were considered.

The proposed algorithm uses SPH, and so the domain is discretized by particles with no connections between them. Through the particles, the approximation functions are derived, and the field function is approximated. Moreover, cell, extracellular matrix, boundary, and blood vessel particles are considered, and the initial velocity, internal pressure, and acceleration of all particles are calculated once the process of proliferation starts. Only one cell was initially considered and was allowed to grow and divide, following an exponential growth. Throughout the whole domain, a VEGF gradient was also defined, and firstly, the highest focus of VEGF was imposed near the cell. As the cells proliferate, a vessel is allowed to grow following the highest VEGF concentration.

Different simulations were run to see the viability of this new feature of the algorithm and to ensure that the generated vessel followed an appropriate path. In every simulation, the vessel grew towards the cluster of cells, which was considered a VEGF source. The growth was also regulated by the cell proliferation process to minimize excessive vessel growth. When this feature was validated, different values of concentration and different points of the focus of concentration were analysed to determine their impact on angiogenesis. Again, in all scenarios, the vessel growth followed the highest concentrations of VEGF, which corresponded to the different focus of concentration in each simulation.

The proposed algorithm was capable of combining the processes of cell proliferation and angiogenesis. As mentioned and reported in the literature, this new feature of the algorithm is dependent on the location of the source of VEGF, which is the factor that triggers the process, and on the cell proliferation process. The vessel growth followed the focus of concentration, as desired and anticipated. Despite the initial phase of the algorithm, the obtained results were reasonable and adequate.

Acknowledgements
The authors acknowledge the funding provided by Ministério da Ciência, Tecnologia e Ensino Superior – Fundação para a Ciência e a Tecnologia (Portugal), under the grant: SFRH/BD/146272/2019, and by LAETA, under the project UIDB/50022/2020.

Neuroblastoma (NB) is the most frequent solid cancer of early childhood. It is a type of cancer that is highly representative of the cancer disease itself, since NB is strongly heterogeneous with very diverse clinical courses that may vary from an indolent disease causing little or no harm and exhibiting spontaneous regression, to an aggressive disease with fatal progression. For these reasons, NB is considered a paradigm of cancer disease and an excellent context of application for the validation of novel developments which have the ambition to be of potential application in a large variety of solid cancers.

NB tumours consist of two main cell populations, neuroblasts and Schwann cells, and the current neuroblastoma classification is based on histological criteria, e. g. the quantity of Schwannian stroma. Neuroblasts and Schwann cells are primary interest herein for contribute directly to the mechanical properties of the tissue through the proliferation and death processes. Extracellular matrix also have a principal role in the cell-microenvironmental cross-talk therefore the tumour can promote to a better stage or keep growing.

We here present a phenomenological model which takes into account as detail as possible to better mimic the real tumour behaviour. Our hypothesis proposes that tumour evolution can be attributed to three distinct processes: growth, shrinkage, and remodelling. The biomechanical model is based on the mass and cellular balance equations coupled with elasticity. The multispecies model simulates the effect of the cellular processes that occur during tumour growth and shrinkage, namely proliferation and death.

The biomechanical finite element model of NB tumour growth starts from imaging data derived mainly from MRI sequences. This data comprises the geometry, the initial cellularity distribution and the tumour vasculature evaluation. At the end of the simulation, the results obtained are validated with a second set of imaging data obtained after treatment.

The study simulates three-month chemotherapy using real patient cases, and presents two distinct outcomes: in one of them, the tumour volume was reduced 20% and in the other one, the volume decreased 90%. One of the patients was classified as low-risk, following the International Neuroblastoma Risk Group (INRG) system, whereas the other was classified as intermediate-risk.
Differences appeared in the histology analysis, which reveal one tumour with a higher concentration of tumoural cells, and in the radiomic data obtained after image analysis. The model effectively reproduces these varying outcomes following the application of chemotherapy, facilitating the identification of cases in which the treatment may be effective.

Acknowledgments

The authors were supported by PRIMAGE (PRedictive In-silico Multiscale Analytics to support cancer personalized diaGnosis and prognosis, empowered by imaging biomarkers), a Horizon 2020 RIA project (Topic SC1-DTH-07-2018), grant agreement no: 826494.
SHR would like to thank the support of the Government of Aragon (Grant no 2019-23).
This work was supported by Grant PID2021-124271OB-I00 and PID2021-122409OB-C21 founded by MCIN/AEI/10.13039/501100011033 and ERDF A way of making Europe.

The Finite Element Method (FEM) has been used in a number of different areas of research, such as orthopedic implant design. Numerical analyses of hip implants and their surface topographies provide important information that can indicate if the connection between bone and implant is good. Although this approach requires less time compared to traditional in vivo studies, it is still time-consuming in situations when a large number of different models need to be created and analyzed. This is the situation with the analysis of a variety of surface topographies, which is needed to better understand the influence of each parameter on shear stress distribution. It is known that higher shear stress values lead to loosening of the contact between femoral bone and hip implant and additional surgical procedures. In order to perform shear stress classification, four algorithms have been considered. Those algorithms were: Support Vector Machines (SVM), K - Nearest Neighbor (KNN), Decision Tree (DT), and Random Forest (RF). A total of 10 model parameters have been used with the previously mentioned classification algorithms in order to obtain information if the shear stress value for the implant will be below the user-defined threshold value (0 – above threshold; 1 – below threshold). The considered parameters were: Number of different radius values (1 or 2); Radius 1 value (>0); Radius 2 value (≥0); Number of half-cylinders lengthwise (>1); Distance between half-cylinders lengthwise (≥0); Number of half-cylinder rows (≥0); Distance between half cylinders rows (≥0); Distance of the first half-cylinder from the edge where loading is defined (≥0); Distance of the last half-cylinder from the edge without loading (≥0); Half cylinder orientation (0 – half-.cylinder follows model length; 1 – half-cylinder follows model width). The database consisted of 60 models, which is the main limitation of this study. The obtained results show that classification algorithms can be useful as a way to have preliminary indications of models that should be further analyzed with FEM. SVM and RF have shown the best results out of the four considered algorithms. For those two algorithms, the following results were obtained: SVM: Precision - 93%; Recall - 92%; F1 score – 92%; RF: Precision - 90%; Recall - 86%; F1 score – 86%. KNN and DT have obtained significantly lower results.

Acknowledgment

This research is supported by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 952603 - SGABU. This article reflects only the author’s The Commission is not responsible for any use that may be made of the information it contains. The research was funded by the Ministry of Education, Science and Technological Development of the Republic of Serbia, contract number [451-03-68/2022-14/200107 (Faculty of Engineering, University of Kragujevac)].

Digital image processing has become increasingly widespread in biological and medical research as microscopy and screening images have become more complex. Automatic image processing plays an important role in increasing the comprehensibility of imaging data inference and reducing the amount of time and manual labor required to fulfill such tasks. One of such tasks requiring digital image processing is Traction Force Microscopy (TFM). Cells adhere to the extracellular matrix (ECM) and exert traction forces to probe the environment, migrate, maintain tissue integrity and form more complex multicellular structures. TFM has become the preferred methodology to quantify exerted forces at the cell–matrix interface and thereby provide quantitative information on cellular mechanics. Before TFM can be utilized, image processing is necessary in multiple phases of data preparation. Data used in this research was acquired from an in vitro model of sprouting angiogenesis in which Human Umbilical Vein Endothelial Cells (HUVECs) invade a polyethylene glycol (PEG) hydrogel. The data is represented as raw 3-dimensional volumes of cells before and after inhibiting their mechanical activity by means of Cytochalasin D. Currently the image processing for TFM required time consuming manual tuning of image denoising filters and tuning of an image segmentation threshold. The dataset was composed of 125 3D volumes of varying sizes, in height and width as well as in depth. These volumes consisted of 9126 total 2D slices, of which 20% was used as the test set, while 80% was used for neural network training.

In this paper we propose a methodology for automatic image segmentation which alleviates most of the manual labor required from the user. The aforementioned stacks of images contained noise caused by either the equipment used or the environment. This noise had to be reduced for the purpose of increasing the quality of training images and segmentation accuracy. The individual 2D slices of each image stack were introduced into the convolutional neural network, while the previously manually segmented images were used as training masks. The segmentation model is based on the U-Net architecture which was modified for this purpose. The model provides 2D segmented images which are concatenated into a 3D stack to reconstruct the 3D input volume. The model achieved admirable results showed by the intersection over union metric, which is greater than 80%. Future research will focus on the improvement of segmentation results through the creation of new image filtering mechanisms and the introduction of more data.

Acknowledgements: The research was funded by the Ministry of Science, Technological Development and Innovation of the Republic of Serbia, contract number [451-03-47/2023-01/200378 (Institute for Information Technologies, University of Kragujevac)]. This research is also supported by the project that has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 952603 (SGABU project). J.B-F. acknowledges the support from the Research Foundation Flanders (FW) grant n. 1259223N. This article reflects only the author's view. The Commission is not responsible for any use that may be made of the information it contains.

In the era of personalized medicine, early and accurate prediction of individuals at high risk of severe consequences of carotid artery stenosis (CAS) which is caused by the atherosclerotic plaque deposition and progression, would allow preventive, therapeutic, or surgical measures before any of life-threatening events take place. It is therefore important to integrate machine learning techniques and computational modeling that can automatically and objectively segment and afterwards mesh the carotid arteries from image data. Among various diagnostic techniques, the US is usually initially recommended CAS diagnostic examination. The US images are used in this study as this diagnostic procedure is low-cost and present at the first level of clinical care, while its improvement in view of image processing and extracted features can be beneficial for more detailed and faster analysis of patients.

Comparing with traditional Computer-Aided Diagnosis (CAD) systems where feature selection and extraction are important steps, deep learning replaces generic imaging features with processing layers that are more complex and more specific to the data. In this study, we have investigated a deep learning approach, which does not require intensity thresholds and imaging features, leading to an optimal use of information and improved prediction power. The automatic carotid artery segmentation was done using U-Net based deep convolutional network (CNN). Deep learning is combined with meshing techniques to perform 3D reconstruction of patient-specific carotid bifurcation, including extraction of atherosclerotic plaque within the arterial wall. The clinical dataset of US anonymized images was used to train and validate the U-Net CNN in automatic segmentation of carotid components (lumen, arterial wall, atherosclerotic plaque) taking into account all three arterial branches (common, internal and external carotid artery). Common classification metrics are considered for one-vs-all quantitative evaluation, including precision (P), recall (R) and F1-score. Results on the test dataset for lumen: P= 0.89, R=0.78, F1=0.83; and for wall: P=0.82 R=0.91, F1=0.85. The results of segmentation are coupled with the reconstruction software to perform a full 3D reconstruction of carotid bifurcation and generate patient-specific 3D mesh. Using the reconstruction software, it is possible to obtain 3D volume information and perform the visualization of carotid structures in different planes and from different angles.

In summary, the automatic extraction of US features related to atherosclerotic carotid artery gives the specific segmentation of individual patient-specific anatomy and can be efficiently integrated with 3D reconstruction. The obtained 3D models can be used for further computational simulations and analyses of individual patients. Integration of deep learning techniques and computer-based modelling can contribute to better risk stratification and assessment of asymptomatic patients with carotid atherosclerotic disease.

Acknowledgements: This paper is supported by the European Union’s Horizon 2020 research and innovation programme TAXINOMISIS (Grant Agreement 755320). This article reflects only the author's view. The Commission is not responsible for any use that may be made of the information it contains. The research is also funded by Serbian Ministry of Education, Science, and Technological Development, grants [451-03-47/2023-01/200378 (Institute for Information Technologies, University of Kragujevac)] and [451-03-47/2023-01/200107 (Faculty of Engineering, University of Kragujevac)].

Mathematical models for multiscale phenomena are typically based on the deduction of a set of differential equations relating the behavior at different involved scales. The number of necessary parameters and the complexity of the equations are crucial in developing effective, predictive models.

A large scientific effort in deducing sophisticated numerical techniques, in particular neural network approaches, characterized the recent literature in the field. This is due to the availability in different fields and in material science of large new data sets down to the nanoscales. New sophisticated approaches exhibit an incredible ability of fitting experimental data even when strongly multiscale nonlinear behavior is observed. The important drawback of these approaches is that typically they are based on statistical based analysis without explicit causal relations, resulting in a poor physical insight.

In this perspective, we propose the adoption of symbolic data modelling techniques, that generate mathematical expressions to fit set of data points using the evolutionary process of Genetic Programming (GP). It is ``Evolutionary Polynomial Regression’’ (EPR) approach integrating regression capability and GP paradigm.

As a result, we obtain explicit analytic formulas that represent the multiscale phenomenon model, elucidating the role of dependent and independent variables, with the important possibility to minimize the complexity of the model in a classical Pareto front.

Specifically, by focusing on the important case of hierarchical biological materials, in this work we deduce numerical techniques that recognize the role of variables in multiscale modelling. Based on a recent set of data at the micro, meso and macroscale of spider silks of different spider species, by means the EPR approach we obtained a set of equations deducing the dependence of the macroscopic thermo-hygro-mechanical behavior from low scales parameters. This set of equations
represents the multiscale model, obtained from data modelling, that we propose to describe the spider silk behavior. Such material behavior constitutes a prototypical example of a physical phenomenon deeply based on a hierarchical behavior and on a complex energetic exchange among different scales. Spider silk is an incredible example on how, based on very `weak’ material components and chemical bonds at the micro scales, it is possible to obtain material behaviors unattained by artificial materials. This is interesting not only for biological aspects, but also in the field of designing new bioinspired materials.

We propose a preliminary analysis of the possibility of adopting such techniques in material science, that in our opinion can represent an important extension of presently adopted numerical methods. As we show, this can be fundamental to gain new physical and modelling insight based on the availability of large new data sets down to the nanoscales. Moreover, we suggest that, due to the generality of our results, the proposed approach can be of interest for many other different fields where multiscale phenomena underly the observed behavior.

In recent years, the use of Artificial Intelligence (AI) models that are informed by physics has become a popular approach for creating data-efficient models in fields such as computational physics, biomedicine, and biomechanics. At the heart of these applications is the complex and adaptable nature of living matter. The adaptations occurring in living systems, such as inflammation, or mechanical changes like growth and remodeling are governed by complex cell-signaling regulatory networks and occur on multiple spatial and temporal scales.

Mechanical modeling of changes in living systems such as peripheral artery plaque progression involves constructing computational models that simulate the development and progression of plaque in the arteries. While these models offer the potential to improve understanding of the underlying mechanisms of plaque formation and to predict the progression of the disease, there are several bottlenecks that can limit their accuracy and effectiveness.

One of the primary challenges is the limited availability of high-quality clinical data including measurements of arterial geometry, blood flow characteristics, and properties of the arterial wall and plaque, all necessary for constructing high quality models. Obtaining this data is often difficult and time-consuming, thus limiting the accuracy and validity of the resulting models.

Creating models that accurately represent complex interactions such as fluid dynamics, solid mechanics, and cellular processes of biological phenomena can be challenging, and models that are too simplistic may fail to capture key aspects of the disease process. In addition, the behavior of arterial plaque is highly variable between patients, and this makes it difficult to develop models that accurately capture the full range of possible disease trajectories.

Finally, creating and running accurate models requires significant computing resources, which can limit the number of simulations that can be performed and can make it difficult to perform sensitivity analyses and explore the effects of different input parameters.

AI has the potential to significantly enhance the accuracy of peripheral artery plaque progression models. The choice of AI method depends on the specific task, available data, and considerations around patient privacy and data security. Machine learning (ML) techniques, including both supervised and unsupervised learning, can be used to predict plaque progression and identify underlying patterns in the data. Genetic Algorithms optimize model parameters for individual patient characteristics and reinforcement learning can determine optimal treatment plans, learning to maximize the reward of minimizing plaque progression. Natural Language Processing (NLP) can be used to extract relevant clinical information from unstructured data, such as clinical notes. Bayesian Networks infer the likelihood of plaque progression based on various patient features. Artificial Neural Networks (ANNs) can model complex relationships between multiple variables, while convolutional and recurrent neural networks, can analyze images or sequential data related to plaque progression.

Acknowledgements

This paper is supported by the DECODE project (www.decodeitn.eu) that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 956470. This article reflects only the author's view. The Commission is not responsible for any use that may be made of the information it contains.

We simulate biophysical processes at various spatial and temporal scales in order to investigate muscle activity. Multiscale simulations, where Huxley's muscle contraction model defines the microscopic scale material properties of muscle, and continuum muscle mechanics is modeled using the finite element approach, can require a significant amount of computational power. The computations performed at the microscale are the most time-consuming component in these simulations. In our work, we developed a computationally more efficient surrogate model to replace the actual Huxley muscle model in order to reduce the computational demands of the simulations.

Instead of solving Huxley’s muscle contraction equation via a method of characteristics, we train a physics-informed neural network to give an approximate solution. Physics-informed neural networks can produce solutions faster than traditional numerical solvers, although these solutions can be less reliable. Once the neural network is successfully trained we use it in multi-scale simulation instead of the traditional method of characteristics. The neural network takes as an input the position of the nearest available actin-binding site relative to the equilibrium position of the myosin head, time, activation, current and previous stretch and produces the probability of cross-bridge formation. For each of the integration points within the finite element model, a large number of positions of available actin sites are observed and stresses are calculated by summing all of the calculated probabilities of the cross-bridge formation.

In our work, we also present the interface between finite elements, at the macro level, and physics-informed neural network at the micro level. The neural network architecture and weights are loaded at the beginning of the finite element simulation, along with the initialization of the neural network input tensor. The size of the neural network input tensor is dependent on the number of integration points in the model and division along the axis at which the positions of the nearest available acting sites are observed. Values of the observed positions are constant during finite element simulation, so only the values for time, activation, and stretches are changed inside the neural network input tensor. To lower time consumption, values for all of the integration points are filled-in, before the neural network’s prediction. Once the tensors are filled-in, the neural network predicts output values for all the integration points in the model. After the neural network has made its prediction, stresses are calculated based on the probabilities of the attached myosin heads, and the finite element procedure continues with the next iteration or next time step. Our physics-informed surrogate model is around two times faster than the original Huxley model solved by the method of characteristics.

This research was supported by the European Union’s Horizon 2020 research and innovation programme under grant agreement No 952603 (http://sgabu.eu/). This article reflects only the author's view. The Commission is not responsible for any use that may be made of the information it contains. Research was also supported by the SILICOFCM project that has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 777204. This article reflects only the authors’ views. The European Commission is not responsible for any use that may be made of the information the article contains. The research was also funded by the Ministry of Education, Science and Technological Development of the Republic of Serbia, contract numbers [451-03-68/2022-14/200107 (Faculty of Engineering, University of Kragujevac) and 451-03-68/2022-14/200378 (Institute for Information Technologies Kragujevac, University of Kragujevac)].

Restenosis refers to the uncontrolled growth of tissue in vessel walls as part of an inflammatory response that follows cardiovascular interventional procedures including balloon angioplasty and stent implantation. Although the risk of restenosis has reduced with the advent of drug-eluting stents, it is not completely eliminated. An in silico replication of neointimal hyperplasia, the mechanism behind restenosis, shall therefore provide the necessary means to derive insights about the biochemical and cellular interactions within the vessel wall, and eventually address the risks of restenosis in a patient-specific manner. In this regard, six important biochemical species in the vessel wall are modeled within the scope of this work: platelet-derived growth factor (PDGF), transforming growth factor (TGF)-β, extracellular matrix (ECM), endothelial cells (ECs), drug (rapamycin) and smooth muscle cells (SMCs). This is considered an extension of the work presented by the authors in [1]. The intricate interactions between these species are then coupled to a continuum mechanical finite growth theory, where the local SMC density drives the growth process and the local ECM (hence collagen) concentration controls the compliance of the vessel wall. Multiphysics-based frameworks for modeling damage-driven growth and remodeling have been presented in several earlier works (e.g. [2]), but the presented model specifically addresses the inflammatory response due to endothelium denudation in addition to the pharmacological influences of the rapamycin-based drugs embedded in modern drug-eluting stents.

[1] Manjunatha, K., Behr, M., Vogt, F., Reese, S., “A multiphysics modeling approach for in-stent restenosis: Theoretical aspects and finite element implementation”, Computers in Biology and Medicine, 150:106166, (2022).

[2] Escuer, J., Martinez, M. A., McGinty, S., and Pena, E., “Mathematical modelling of the restenosis process after stent implantation”, J. R. Soc. Interface, 16:20190313, (2019).

Aortic Dissection is a severe medical condition that can lead to fatality if not treated early. Aortic Dissection begins when a tear occurs in the aortic wall's inner layer (intima). Moreover, a second blood-filled channel called a false lumen is created where thrombosis most probably will happen. In type B Aortic Dissection, severe complications, such as malperfusion or rupture in the aorta, can occur in the first seven days, resulting in high mortality rates. The optimal treatment for type B Aortic Dissection remains a topic of debate, leading to an interest in computational methods to aid in decision-making for treatment.
To obtain a macroscopic model the Theory of Porous Media is chosen to modell the multiphasic structure of the thrombus. The whole aggregate is divided into solid, liquid, and nutrient constituents. As usual in the theory of porous media, all three phases occupy each material point simultaneously. The constituents are assumed to be materially incompressible, under isothermal conditions, and the entire aggregate is fully saturated. Darcy's law is used to describe the flow of fluid phases, i.e., the liquid and the nutrient phase. The constituents' balance equations have coupling terms that account for the interactions between the different phases. The essential coupling term for the growth process is the mass exchange term in the volume balances. Due to the incompressibility assumption the volume fractions can be used to define the constituents. Therefore, the regions with thrombus are determined using the solid volume fraction.
Initially, a set of equations for thrombosis in type B Aortic Dissection is presented. The respective variational form is approximated with finite elements and implemented in PANDAS, a finite element package used for solving strongly coupled multiphase porous media problems. However, the model's description is challenging due to the need for detailed knowledge and parameters to quantify the influence of different factors. Nevertheless, the effects of blood velocity and nutrients on thrombus growth are well researched. Therefore, a velocity and nutrient concentration induced growth model is proposed within the Theory of Porous Media. The effects of various parameters are studied using a realistic 2D geometry of the false lumen, and reasonable parameters for thrombus growth are selected. We present the model's features in agreement with the Virchow triad and its usefulness in actual cases. The simulations show that the thrombus grows in the low-velocity regions of the blood. The proposed model provides a reasonable approach for the numerical simulation of thrombosis.

Endovascular coil embolization is one of the primary methods for the treatment of a ruptured or unruptured cerebral aneurysm. Though being a well accepted and minimally invasive method, it bears the risk of suboptimal coil placement, which can lead to incomplete occlusion of the aneurysm, causing a degenerate healing process or in the worst case even aneurysm recurrence.

Our study investigates factors that lead to a suboptimal placement when a coil is deployed in an aneurysm. To this end, we use a mathematical model that simulates the embolization process of the coil in patient specific three-dimensional aneurysm geometries. One of the key features of coils is that they have an imprinted rest shape, which helps to safely fix them inside aneurysms. We consider the latter property by introducing bending and torsional forces into our model of the coil. In order to differentiate optimal from suboptimal placements, we use measures like the local packing density, which is connected to the quality of the subsequent healing process.

We model the coil as a one dimensional discrete space curve with the dynamics of a mass spring system. For this purpose, we include stretching, bending and torsional springs. The aimed for deflected rest shape is accounted for by penalizing deviations from it within the individual spring energies. To obtain the forces that act on the coil, we differentiate the total elastic energy of the spring mass system by algorithmic differentiation with respect to the material points. Collisions between coil segments and the aneurysm-wall are taken into account by an efficient contact algorithm, such that only a fraction of the possible contact partners have to be evaluated. The numerical solution of the model is obtained by a first order accurate time stepping scheme.

We then carry out a sensitivity study that shows how process parameters like the insertion point or the angle of insertion affect the shape of the deployed coil in the aneurysm. We further show how other medical tools such as, for example, stents or hyper-elastic balloons can support the coil-placement process. Finally, we investigate the sensitivity of our model with respect to the material parameters and rest shape of the coil.

Our model can be easily incorporated into blood-flow simulations of embolized aneurysms. It is our goal to eventually provide a simulation tool that will provide new insight into the behavior of coils and allow medical professionals to better plan coil embolization procedures and hence assist them in their practice.

Left ventricular (LV) hemodynamics are hypothesized to serve as an early indicator for the manifestation of cardiovascular diseases. In the analysis of hemodynamics, computational fluid dynamics (CFD) can complement medical imaging methods for in-depth flow investigations. To obtain information about individual conditions, a high level of patient-specific input into the model is required. Considering the variety of data sources, medical imaging modalities, and the range of data quality, this leads to demanding case processing with a lot of manual effort.

To target this issue, we standardized the workflow by coupling the patient-individual anatomical representation in form of a statistical shape model (SSM) with a fluid structure interaction (FSI) workflow. Thereby, we obtain a high degree of individuality in our simulations while keeping the case-specific effort involved at an acceptable level. Furthermore, the workflow allows usage of all common cardiac imaging modalities, namely computed tomography (CT), magnetic resonance imaging (MRI), and ultrasound (US).
To test the simulation pipeline, we processed data of healthy volunteers by comparing the FSI-computed flow results to measured velocity fields from 4D flow MRI. Further applications of the workflow were tested with simulations of patients showing isolated or combined aortic stenosis and mitral regurgitation. For patients with ventricle aneurysms after myocardial infarction, we compared pre- and post-treatment situations.

Imaging data of the healthy volunteers included 4D flow MRI and cine MRI sequences. While the latter served as model input, the flow data was used to check the results of computed hemodynamics for plausibility. Simulations of mitral regurgitation patients revealed characteristics of altered kinetic energy and specific energy dissipation during diastole. These results correspond well with findings of previous clinical MRI studies. Simulations of the pre- and post-surgical hemodynamics of left ventricular aneurysms showed significantly improved flow conditions after surgery. Although, the flow results cannot quantitatively be compared to clinical data because no MRI was acquired, results show improved flow patterns and align well with the improvement of the medical conditions patients showed after surgery.

To conclude, we could prove that our workflow is robust and easy to use while including a lot of patient-specific details. This allows for representation of individual flow features, which can perspectively be used in the clinical routine to support diagnostics, therapy planning, and outcome prediction. A profound validation study must be performed to test the approach on a sufficient amount of data. However, the availability of such comprehensive data is rare and prospective studies are very costly.

The left atrial appendage (LAA) is a small, finger-like structure in the heart's left atrium. Despite its small size, the LAA has been shown to affect cardiovascular health significantly. One of the challenges associated with the LAA is that it is its location, an area of the heart with poor blood flow, which can lead to blood stagnation. Patients with atrial fibrillation often experience insufficient contraction of the left atrium, predisposing the LAA morphology to hemostasis and thrombus formation, resulting in an increased risk of cardioembolic events. To prevent these pathologies, oral anticoagulation therapy is typically used as the primary treatment option for patients. However, not all patients are eligible for long-term oral anticoagulation therapy, which can cause complications such as bleeding. These circumstances led to the development of alternative treatment options, such as percutaneous occlusion devices. However, several drawbacks remain. Peri-implant leakage and device-related thrombosis are common complications in LAA closure procedures. Efforts have been made to reduce these risks, but interpatient heterogeneity remains challenging. Ongoing research aims to develop better treatment options for patients with atrial fibrillation and other cardiovascular conditions. One area of innovation is additive manufacturing (AM), also known as 3D printing, which can improve the accuracy of the selection of LAA closure devices. AM allows for the creation of complex and precise structures with high levels of customization, making it an attractive tool for personalized medicine. AM can generate personalized LAA used in medical practice device simulations, reducing the risk of device-anatomy mismatch and improving the procedure's success rate. Additionally, numerical simulation techniques are being employed to model the behavior of LAA, allowing researchers to optimize device performance and minimize the risk of complications.

This study presents a structural numerical approach for analyzing the optimal material for 3D printing a personalized LAA through finite element analysis (FEA). A 3D model of an LAA was obtained from an actual patient's computerized tomography (CT) scan and subsequently modeled numerically. The material selection for computational analysis is crucial to mimic human tissue's mechanical behavior accurately. When designing a cardiovascular training model, the chosen material for the LAA model must withstand the LAA pressures of device implantation. Higher radial resistance often correlates with higher tensile resistance in materials. Mechanical tensile tests were performed to evaluate the radial resistance of multiple materials. Thermoplastic polyurethane (TPU) material exhibited significant deformation during the tests, reaching 40% without breaking, which indicates its potential as a suitable material for replicating the biological tissue of the LAA. The ADINA® software was employed to perform the finite element analysis (FEA). The study yielded a maximum displacement of 0.4 millimeters for the LAA model, demonstrating a close resemblance to a real LAA in FA condition. Therefore, it was possible to conclude that TPU is a potential material to produce the LAA model for pre-procedural occlusion planning. By reducing the risk of complications associated with percutaneous occlusion devices, these innovations can improve the outcomes and quality of life for patients with AF and other cardiovascular conditions.

The range of materials available for 3D printing has been expanding rapidly in recent years. However, for very specific requirements, such as an anatomical model, there is not always a suitable material available. In particular, the design of gradients in elasticity, color, or surface properties is not truly represented by pure materials. Material jetting allows 3D printing of multiple materials simultaneously, resulting in composite materials with new properties. This study investigated and compared the mechanical properties of pure and composite materials and the possibilities of predicting composite properties by knowing the proportions of the pure materials used.

Samples of commercially available materials (VeroClear, RGD8530, Stratasys Ltd., Minnesota, USA) in their pure and mixed matrix inclusion forms were produced using material jetting (Connex3, Stratasys Ltd., Minnesota, USA). The composites had 1mm³ unit cells, including a cubic inclusion with a volume fraction of f_inc=10%, 30% or 45% RGD8530 in a VeroClear matrix. They were mechanically characterized by uniaxial tensile testing according to ISO 527 for Young’s modulus and Poisson’s ratio. In order to find and validate a method for predicting the properties of the multimaterials, multimaterial homogenization and finite element (FE) modelling were evaluated and compared with the measurement results. Inclusion size and geometry were characterized by optical coherence tomography (OCT) and digital microscopy.

The materials had Young's moduli ranging from 800MPa to 2.5GPa. Multimaterial composites were never as stiff as the primary materials' volume weighted average (26.5±2.7% softer for 45% inclusion volume). OCT scans revealed deviations from the digital design, more specifically a rounding of the inclusion edge, as well as blurred material interfaces within the polyjetted layers.

Models that assume ideal interface conditions, such as multimaterial homogenization or conventional FE simulation, are not capable of producing the measured Young’s moduli. A functional simulation model to predict the uniaxial Young’s modulus, can be established using FE simulation and considering a contact stiffness of F_A=2.2TN/m³ and an inclusion edge radius of r=220, as seen in the OCT data.

In conclusion, matrix-inclusion composites have a non-trivial behavior of elastic properties, which needs an adapted model to predict. The established FE model incorporates the interface stiffness between the individual materials used and the geometric deviations from the digital design that occur during the 3D printing process. This enables more complex parts to be produced using less primary material by predicting, designing and 3D printing structures with precisely defined mechanical properties and gradients that are required to mimic biological structures.

Background

Three-dimensional (3D) digital and additively manufactured models are increasingly used for pre-operative planning, especially for orthopaedics applications. However, the size of the model’s minimal detectable features in clinical CT imaging has not been determined, hence it remains unknown which features might remain undetected. Furthermore, it is crucial to identify potential error sources during CT imaging and image processing, quantify their effect on the 3D model accuracy, and develop an optimized workflow to allow for the reliable use of 3D models in the clinical routine. This study aimed to investigate the minimal detectable bone feature size in CT images and corresponding digital 3D models and to determine the errors attributed to different CT technologies, scanners, scan protocols (clinical versus high dosage for improved model quality), segmentation algorithms, as well as specimen orientation and consequently quantify the resulting geometrical deviations on a defined bone fracture model.

Materials and Methods

Incisions in the diaphyseal radii with 200 and 400 μm width, and bone lamellae (bony displacements) with 100, 200, 300, and 400 μm width were generated in twenty paired forearm anatomic specimens (age: 78 ± 8 years (5 male and female, each)). Additionally, a throughout osteotomy was performed in the diaphysis, held in place with additively manufactured guides to simulate a complete 100 μm wide fracture. Specimens were scanned with different CT scanners and corresponding digital 3D models were created. The effects of CT technologies/scanners, specimen positionings, scan and segmentation protocols, and image post-processing settings on feature detectability were assessed. Furthermore, the intra-and inter-operator variabilities were assessed for the segmentation threshold and 3D model accuracy. Three-dimensional reconstructions of surface scans of the physical specimens were used as ground truth to compute the specific geometry deviation.

Results

In CT images, fracture gaps of 100 μm, and bone lamellae of 300 μm and 400 μm were identified at a rate of 80 to 100%, respectively, independent of the investigated settings. In contrast, only 400 μm incisions and bony displacements were visible in digital 3D CT models. Hereby, the detection rate was independent of the scan settings but dependent on the selected CT technology. image segmentation and post-processing algorithms. Intra- and inter-operator variability were fair to excellent for 3D model accuracy (intra-class correlation of 0.43 to 0.92) and mean 3D deviation was < 0.16 mm on average for all operators, using a simple global segmentation threshold and minimal post-processing.

Conclusion

This first systematic investigation of the effect of multiple variables affecting 3D model accuracy demonstrated that sub-voxel imaging resolution was achieved for all variables. Thus, state-of-the-art CT imaging allows for the detection of bone features down to 100 μm. Corresponding digital 3D models still enable the identification of 400 μm features but require verification with the original CT image series.

Acknowledgements

This work has been partially supported by the Austrian Research Promotion Agency with the project “Additive Manufacturing for Medical Research, M3dRES (nr: 858060).

Hallux Valgus, also called a bunion, is a deformation of the metatarsophalangeal joint, more common in adults, with an estimated prevalence of 23% in the 18-65 age group and 35.7% over 65 years. There is a 3:1 relation for females, tending to increase with age [1]. Several reasons contribute to this type of deformity associated with tight shoes and high heels [2]. One of the theories is based on the idea that a continuous overload leads to deformation of the first metatarsal, with stretching of the capsule. The resulting imbalance of the first metatarsal phalangeal joint leads to increased deformity. The condition worsens, forming a bulge on the medial side [3].

The corrective treatment of the Hallux Valgus deformity usually involves surgery, with several surgical procedures reported. However, the Chevron Osteotomy is the procedure usually followed in this type of treatment, consisting of translating a distal portion containing the head of the first metatarsal, produced by the cut in the sagittal plane of a 60º angle of the distal apex. The resulting central exostosis is excised, with subsequent plication of the medial capsule [4]. This surgical procedure is performed manually, without support systems, depending on the result of the surgeon's experience.

This work proposes a new fixation device that allows the blade's positioning, stabilization, and guidance to guarantee precision in the cutting procedures in the Chevron Osteotomy. The development of the device involved the identification of the criteria for positioning the blade and necessary positional adjustments to ensure the best alignment and stabilization of the cutting blade to perform the Osteotomy. The methodology starts with the geometry of foot bones, obtained from a computed tomography scan using the Mimics Innovation Suite^Ò software. Based on the geometry, a modular support device was conceived and 3D modelled (with Solidworks^Ò software) to be anchored in two Kirschen wires to be applied, one in the first metatarsal and the other in the first proximal phalanx.

The device was prototyped, and a protocol for the experimental tests was defined. Artificial bones were produced from the geometry of the bones, considering the cortical and trabecular components. The tests were conducted in the laboratory with an experimental simulation of the Chevron Osteotomy, varying the positional adjustments associated with the technique. The obtained results demonstrate the device's effectiveness with well-defined cuts.

References

[1] Nix S., Smith M., & Vicenzino B. Prevalence of hallux valgus in the general population: a systematic review and meta-analysis. Journal of foot and ankle research, 3, 21. https://doi.org/10.1186/1757-1146-3-21, 2010.

[2] Sharma J, Arora A, Gupta S. Disorders of the toe - Hallux valgus. Textbook of Orthopaedics & Trauma, Vol. 4. Jaypee Brothers Publishers, 2019.

[3] Lerat J-L. Cheville-Pied: L’Hallux Valgus. Traumatologie Orthopedie. Centre Hospitalier Lyon-Sud, Service de Chirurgie Orthopedique et de Medecine du sport, Cap. 6, 2011.

[4] McKean J., Park J. Hallux Valgus. Lineage Medical, Inc. https://www.orthobullets.com/foot-and-ankle/7008/hallux- valgus, janeiro 2023.

Acoustic waves propagation in biological tissues has been studied usually in the context of diagnostic methods, such as the elastography. In microfluidic devices, nonlinear acoustic phenomena, namely the acoustic radiation, acoustic streaming are employed to manipulate particles and actuate fluid flow. These principles have attracted much interest of the research focused on developing new tissue engineering technologies since the acoustic wave are highly biocompatible, providing a non-contact controlable handle to manipulate bioparticles, or cells.
The paper reports on development of multiscale models of the acoustic waves and related phenomena in porous fluid saturated structures representing perfused tissues, or biomaterial structures intended for fluid suspensions and particle transport due to acoustic power, taking into account nonlinear interactions between to deforming pores and particles.
We focus on modelling two important phenomena induced in porous structures by acoustic waves: the peristaltic deformation and the acoustic streaming. Both can propel the fluid and particle transport, influencing interaction between particles and pore walls. This problem is motivated by the recellularization in decellularized tissue scaffolds, an important progressive tissue engineering technology .
The fluid-structure interaction problem is imposed in elastic scaffolds. To capture both the peristaltic deformation and acoustic streaming in response to propagating acoustic wave, nonlinearities originating in the divergence of the Reynolds stress, the advection acceleration term in the Navier Stokes equation, and the nonlinearity generated deforming pore geometry must be retained. The perturbation with respect to a small parameter proportional to the inverse Strouhal number is applied. This yields the first and the second order sub-problem enabling to linearize the Navier-Stokes equations governing the barotropic viscous fluid dynamics in periodic scaffolds. Subsequent treatment by the asymptotic homogenization leads to a two scale problem where the macroscopic model of the porous medium describes the acoustic streaming. The vibro-acoustic analysis of the first order problem yields the streaming source term for the second order problem which attains the form of the Biot-type medium.
To respect the peristaltic deformation wave in the homogenized model derived in the fixed reference configuration, influence of the deforming channel geometry on the homogenized coefficients must be considered. We use the linear expansions of the permeability and other poroelastic coefficients of the homogenized model; these are based on the sensitivity analysis of the homogenized coefficients with respect to macroscopic strain and pressure which determine the deformation of local microconfigurations.
To account for the acoustic waves influence on the propulsion of cells in the tissue scaffolds, we consider microstructures with trains of inclusions flowing in pores. Using simplified assumption and additional kinematic constraints, the contact of cells with pore walls is avoided. The two-scale homogenization-based modelling provides analytical and computational tools to simulate the influence of the bulk acoustic waves on the microflow of the biological fluids in tissues. In particular, the presented research is aimed at providing a computational feedback to experimental studies of the acoustic waves impact on the tissue scaffolds recellularization. The models are implemented in our in-house developed finite element based software SfePy. Numerical illustrations are presented.

Many liver diseases display an increasing incidence, whereby hepatocellular carcinoma rank 4^th in mortality among cancers. A deep understanding of liver functioning and disease requires to study the complex orchestration of processes at many levels of organization, which favor a systems medicine approach, where experiments, clinical data acquisition and computational modeling are integrated.

The computational models are able to verify or falsify hypothesized mechanisms and already led to the identification of unrecognized mechanisms that could be used to propose new therapeutic strategies [Ghallab et. al., J. Hepat. 2016]. Along the same line of argument, a systems medicine approach recently questioned the consensus mechanism of bile transport by flow in bile canaliculi, the smallest bile conduits in the liver [Vartak et. al. Hepatology 2020], and demonstrated that novel experimental results were compatible with diffusive transport alone. The latter work focused on transport in the bile canaliculi network but did not take into account the entire pathway of bile salts from the liver entrance into the liver capillaries (sinusoids), through the hepatocytes, their excretion into the bile canaliculi and then their flux to the bile ducts from where they are transported to the gall bladder.

Here we will present the necessary modeling components for such a model at the level of liver tissue micro-architecture, from which some have recently been established, as, e.g., blood hemodynamics in sinusoids [Boissier et. al. 2021], and flux of metabolites into and out of hepatocytes [Dichamp et. al., 2023]. The works vary a number of parameters, e.g., in [Boissier et. al.] different numerical algorithms are compared with regard to accuracy and computational efficiency. For clinical relevance, such an entire circulation model shall study the effect of heterogeneity, and permit scalability.

Given the complexity, these requirements may make such models computationally inefficient or even infeasible. To overcome this, various model reduction techniques can be applied, that partially depend on the nature of available data.

We present recent research and first steps towards a proof-of concept model of fluorescent marker transport in the blood-hepatocyte-bile system at the level of liver micro-architecture displaying each individual cell. Starting from general 3D model of mixture flow and transport in each part of the system, we present a strategy for reducing the governing equations into 1D. We then follow by formulating a mathematical description of exchanges occurring at compartment interfaces. Numerical solutions of the model are presented for cases corresponding to used experimental setups. Even with the rather oversimplified assumption on the system geometry (1D) in the prototype model, we are able to capture qualitatively the zonation features that are observed experimentally in mice.

In the future, this model will be refined to explain image data on vascular and biliary networks obtained in intravital imaging. An extension towards upscaling could start from interpreting the parameters of the model at cellular resolution as effective values of their real-world counterparts and go to more complex techniques, such as homogenization (cf. Rohan et al. [2021a], Rohan et al. [2021b]).

The flexoelectricity corresponds to the couplings between the electric polarization of a solid and strain gradient. Identified in the sixties, this nanoscopic electromechanical interaction presents two main attractive features. Firstly, this phenomenon may exist in materials with centrosymmetry, contrary to the piezoelectric effect. Secondly, strain gradient effects being inversely proportional to the medium size, the scale of flexoelectricity is then rather large. Flexoelectricity still presents challenging experimental and theoretical issues [1] and offers promising avenues of research and promises as new therapies in the context of regenerative medicine [2]. Up to now, recent advances in the understanding of bone physiology point out the possible relevance of flexoelectricity effects in bone. Indeed, putting forward the seminal studies of Williams [3], Vasquez-Sancho et al. [4] experimentally showed that part of bone electricity may be due to the bone’s mineral flexoelectric property in addition to already described collagen’s piezoelectricity and stress-generated streaming potentials. This may explain the in vitro observation of osteoblasts (bone forming cells) activity near a crack in pure hydroxyapatite (HAP), without the presence of collagen (COL) [5]. Since this electrical polarization and its consequences in interstitial fluid flow within bone is a key issue in its adaptation and repair [6,7], the high potential of smart flexoelectric tailored biomaterials in the highly competitive osteo-articular tissue engineering domain has to be investigated.

In this context, this study aims at correlating the crack-induced flexoelectric effect combined to collagen piezoelectricity on the bone cell behaviour in its fluid environment. It is thus necessary to connect the microscopic phenomena in the vicinity of bone cells to bone tissue structures. Thus, a trans-scale approach (homogenization) is carried out to propagate nanoscopic electromechanical consequences at the upper scale where the bone modelling units (BMU) activity occurs. In particular, we derive from this multiscale approach the coupled phenomena governing the osteocytic sensing and cell-to-cell communication known to drive BMU. This study potentially may have in a longer term a great technological impact with the design of novel biomaterials and technological solutions to stimulate bone repair.

[1] Chae et al., ACS Appl Bio Mater, 1: 936, 2018.

[2] Wang et al., Prog Mat Sci, 106:100570, 2019.

[3] Williams & Breger, J. Biomech, 8:407, 1975.

[4] Vasquez-Sancho et al., Adv Mater, 30:1705316, 2018.

[5] Shu et al., Mater Sci Eng C, 44:191, 2014.

[6] Lemaire et al., J Mech B Biomed Mat, 4:909, 2011.

[7] Lemaire et al., Int J Num Meth Biomed Eng, 29:1223, 2013.

Adherent cells in vivo often reside on basement membranes which are thin, sheet-like specialized extracellular matrices that function as cellular anchorage sites, physical barriers, and signaling hubs. An important physical attribute of basement membranes is their three-dimensional topography due to the complex microstructural organization of their constituent proteins. This topography takes the form of an intricate network of fibers and pores with characteristic dimensions on the order of 100 nm on top of which lies a larger, micron-scale topography in the form of anisotropic undulations. Various types of engineered microstructured substrates have been developed to study the impact of topographical cues on cellular structure and function in vitro. One such system that we have been using is a substrate that consists of arrays of microgrooves that are intended to mimic the anisotropic organization of the basement membrane and on which cells can be directly cultured.

Our studies of the effect of substrate topography on cellular structure and function have focused principally on the vascular endothelium, the monolayer of cells lining the inner surfaces of all blood vessels. In medium and large arteries, chronic endothelial inflammation is a trigger for atherosclerosis, the disease that leads to heart attacks and strokes. Interestingly, atherosclerotic lesions develop preferentially in arterial regions where endothelial cells are cuboidal and randomly oriented, whereas arterial zones that are characterized by highly elongated and aligned endothelial cells remain largely spared from the disease. Therefore, understanding the relationships between endothelial cell shape/alignment and function is of fundamental interest, a question that we are tackling using topographic microgroove substrates.

In this presentation, I will discuss four aspects of endothelial cell responsiveness to microgroove substrates and highlight opportunities for computational modeling in each of these aspects. First, I will show how microgroove substrates can be used to noninvasively control endothelial cell shape and alignment and will describe our understanding of the mechanisms that underlie cell shape regulation by microgrooves. Second, I will describe dynamic live-cell recordings that demonstrate that microgrooves can orient the direction of migration of endothelial cells within monolayers and can lead to a unique pattern of collective cell migration that takes the form of antiparallel streams. Modeling the endothelial monolayer as an active fluid with the effect of the microgrooves considered as an energetic constraint on cell orientation predicts the emergence of the antiparallel streams as well as the dimensions of these streams. Third, I will show how microgrooves lead to extensive deformation of endothelial cells and their nuclei and will evoke the interesting notion of using these deformations to diagnose certain diseases that involve abnormalities in cellular and nuclear mechanical properties. Finally, I will describe the competition between microgroove-derived contact stresses on the cells’ basal surface with flow forces on the cells’ apical surface and how this competition is a key determinant of endothelial cell shape and alignment.

Osteoporosis (OP) is a chronic progressive bone disease which affects a large portion of the elderly population worldwide. OP is characterized by a slow reduction of bone matrix and changes in the bone matrix properties which ultimately leads to whole (organ) bone fractures [1].

Novel drug treatments are developed to more effectively reduce the risk of bone fractures. Assessing the effects of novel and existing treatments on OP is challenging due to the complexity of the bone remodeling process, its effects on the bone matrix and the different spatial and temporal scales involved. Identification and characterization of various bone biomarkers has significantly improved our understanding of OP pathophysiology. The bone matrix and its constituents are specific bone biomarkers measured at a particular bone site. On the other hand, biochemical ligands released during bone remodeling and measured in blood or urine are non-specific bone biomarkers. These biomarkers can be used to characterize the underlying bone mechanobiological system and drug treatment effects [1].

Recently, disease system analysis (DSA) has been proposed as a novel approach to quantitatively characterize drug effects on disease progression [1]. DSA integrates physiology, disease progression and drug treatment in a comprehensive mechanism-based modelling framework using a large amount of complementary biomarker data. DSA applied to population based structural models of whole bony organs (e.g. femur and vertebra) can be used to perform in silico trials of drug efficacy for osteoporosis treatment. In this presentation, I will present latest mechanistic pharmacokinetic-pharmacodynamic (PK/PD) models of osteoporosis treatments. Examples of currently used drug interventions including denosumab [2,3] romozosumab [4], and PTH [5] treatments will serve as discussion points on which mechanisms are essential for accurate bone remodeling simulations. Bone matrix mineralization turns out to be an essential model feature that is required to predict BV/TV changes for the case of anti-catabolic drug treatments of OP [3]. Finally, I will elucidate on how PK/PD models, typically applied at the tissue scale, can be adapted to the whole organ scale and coupled with structural finite element simulations of bone in order to predict effects of drug treatments on bone strength and fracture risk.

Acknowledgments: Dr Pivonka acknowledges support from the Australian Research Council (IC190100020) and (DP230101404).

References: ”[1] S. Trichilo and P. Pivonka, Disease systems analysis in osteoporosis and mechanobiology, in Multiscale mechanobiology of bone remodelling and adaptation, Editor P. Pivonka, CISM Courses and Lectures No. 1406, Springer, 2017; [2] S. Scheiner et al.. Mathematical modeling of postmenopausal osteoporosis and its treatment by the anti-catabolic denosumab, Int. Journal for Numerical Methods in Biomedical Engineering, 30(1), pp1-27, 2014; [3] J. Martinez-Reina and P. Pivonka, Effects of long-term treatment of denosumab on bone mineral density: insights from an in-sillico model of bone mineralization, Bone, 125, pp87-95, 2019; [4] M. Martin et al., Assessment of Romozosumab efficacy in the treatment of postmenopausal osteoporosis: results from a mechanistic PK-PD mechanostat model of bone remodeling, Bone, 133, pp1-16, 2020; [5] M. Lavaill et al., Effects of PTH treatment in osteoporosis, BMMB, pp1-16, 2020.”

As the key organ for metabolic processes in the human body, the liver is responsible for essential processes like fat storage or detoxification. Liver diseases can trigger growth processes in the liver, disrupting important hepatic function-perfusion processes[1]. To better understand the interplay between hepatic perfusion, metabolism and tissue in the hierarchically organized liver structure, we developed a multicomponent, poro-elastic multiphasic and multiscale function-perfusion model [2,3], using a multicomponent mixture theory based on the Theory of Porous Media (TPM). The multiscale approach considers the different functional units of the liver, so-called liver lobules, with an anisotropic blood flow via the sinusoids (slender capillaries between periportal field and central vein), and the hepatocytes, where the biochemical metabolic reactions take place. On the lobular scale, we consider a tetra-phasic body, composed of a porous solid structure representing healthy tissue, a liquid phase describing the blood, and two solid phases with the ability of growth and depletion representing the fat tissue and the tumor tissue. The phases consist of a carrier phase, called solvent, and solutes, representing microscopic components, e.g. nutrients, dissolved in the solvent. To describe the influences of the resulting tissue growth, the model is enhanced by a kinematic growth approach using a multiplicative split of the deformation gradient into an elastic and a growth part, dependent on the fat accumulation and tumor development. To describe the metabolic processes as well as the production, utilization and storage of the metabolites on the cellular scale, a bi-scale PDE-ODE approach with embedded coupled ordinary differential equations is used. In order to represent realistic conditions of the liver, experimentally or clinically obtained data such as changes in perfusion, material parameters or tissue morphology and geometry are integrated as initial boundary conditions or used for parametrization and validation [4]. Data integration approaches like machine learning are developed for the identification, processing and integration of data. A workflow is designed that directly prepares the model for clinical application by (semi-)automatically processing the data, considering uncertainties, and reducing computation time.

Acknowledgements: This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) via the following projects: 312860381 (Priority Programme SPP 1886, Subproject 12); 390740016 (Germany’s Excellence StrategyEXC 2075/1, Project PN 2-2A); 436883643 (Research Unit Programme FOR5151, Project P7); 465194077 (Priority Programme SPP 2311, Project Sim-LivA); 463296570 (Priority Programme SPP 1158, Antarctica)

REFERENCES
[1] Christ, B., ..., Ricken, T. et al. [2017], Front. Physiol., doi: 10.3389/fphys.2017.00906 .
[2] Ricken, T., et al. [2015], Biomech. Model. Mechanobiol., doi: 10.1007/s10237-014-0619-z .
[3] Ricken, T., and Lambers, L. [2019], GAMM-Mitteilungen, doi: 10.1002/gamm.201900016 .
[4] Seyedpour, S. M, ..., and Ricken, T. [2021], Front. Physiol., doi: 10.3389/fphys.2021.733393

Physiologically-based pharmacokinetic (PBPK) models are powerful tools for studying drug metabolism and its effect on the human body. Here, we present our work on PBPK models for metabolic phenotyping and liver function evaluation [1-5]. To develop and validate our models, we established the first open pharmacokinetics database, PK-DB, which contains curated data from over 600 clinical studies. Our models can be individualized and stratified, allowing us to simulate the effect of lifestyle factors and co-administration on drug metabolism.

We have applied our models to various clinical questions, such as simulating individual outcomes after hepatectomy using an indocyanine green model and studying the effect of CYP2D6 gene variants using a model of dextromethorphan coupled with drug-gene interactions. These models are constructed hierarchically, describing metabolism and other biological processes of organs such as liver and kidney coupled to whole-body physiology. All models and data are available for reuse in a reproducible manner encoded in the Systems Biology Markup Language (SBML).

One of the major challenges in PBPK modeling is the integration of these models with cellular and continuum-biomechanical models. By coupling PBPK models with these models, it becomes possible to provide drug concentrations at the site of action as boundary conditions and time-varying inputs to continuum-biomechanical models. This approach enables us to study the impact of systemic chemotherapy on tumor growth or to investigate changes in perfusion and pressure at the whole-body level on transport and biomechanical properties at the tissue scale.

To demonstrate the effectiveness of this approach, we present an example of PBPK modeling coupled with a multiscale and multiphase continuum-biomechanical model for the liver. Our model integrates partial differential equations (PDE) on the lobular scale with ordinary differential equations (ODE) on the cellular and whole-body scale, resulting in a PDE-ODE coupling. In our showcase, we investigate the drug metabolization process in the liver, considering the possible damage to the liver tissue due to necrosis caused by a toxic side product. The substrates and products of this process can be transported via the systemic circulation and be removed by the kidneys. By coupling PBPK models with continuum-biomechanical models, we can study the transport and biomechanical properties of these substances at both the tissue and whole-body levels. This coupling approach provides a promising method for studying drug metabolism and its effect on the human body at multiple scales, ranging from cellular to whole-body.

Our approach allows us to investigate the impact of this drug metabolization process on the liver tissue and how it affects the overall functioning of the body. Additionally, we can study the impact of various factors, such as different drug dosages, on this process and determine the optimal drug dosage for a given patient. This integrated approach provides a comprehensive understanding of the drug metabolization process and its impact on the human body, paving the way for more effective drug development and personalized medicine.

[1] https://doi.org/10.1093/nar/gkaa990
[2] https://doi.org/10.3389/fphar.2021.752826
[3] https://doi.org/10.3389/fphys.2021.730418
[4] https://doi.org/10.3389/fphys.2021.757293
[5] https://doi.org/10.1101/2022.08.23.504981

Liver transplantation is the only curative treatment option for acute and chronic end-stage liver disease. As a result of demographic change and a more western way of life, the number of elderly multi-morbid prospective recipients and donors grows. Liver grafts from such donors, so-called marginal liver grafts, are often affected by hepatic steatosis compromising the quality of the donor organ. One reason for this is the alteration of the tissue structure, resulting in an impaired perfusion, in turn affecting hepatic metabolism and organ function. In case of a marginal graft, the surgeon is faced with the clinical decision to accept or reject the organ, increasing either the recipient’s postoperative risk or the risk of death on the waiting list significantly. Two major challenges for marginal grafts are the storage between procurement of the organ and transplantation (cold ischemia time) alongside damage after reperfusion, the so-called ischemia reperfusion injury (IRI).

For this purpose, we use computational multiphase and multiscale continuum-biomechanical modeling of biological tissue to simulate the hepatic deformation-perfusion-function relationship (cf. [1, 2]). Based on the theory of porous media (TPM, cf. [3]), we describe the functional liver units, the liver lobules, as a homogenized porous medium. A poroelastic multiphase and multiscale function-perfusion model is obtained by considering the porous liver tissue saturated with an anisotropic blood flow and coupling the metabolic processes at the cellular level within a bi-scale approach (cf. [4]). This approach combines partial differential equations (PDE) on the lobular scale with ordinary differential equations (ODE) on the cellular scale resulting in a PDE-ODE coupling. At the lobular scale, we consider healthy liver, necrotic, and fat tissue as solid phases, while blood is given as a fluid phase. Based on this, systems biology models can be used to model the energy balance, cell death, and functionality of each hepatocyte at the cellular level. Thus, it is possible to describe the ischemic damage caused by nutrient depletion, in turn influencing the perfusion at the lobular scale.

To make patient-specific predictions, we enhance the model through the integration of experimental and clinical data for validation and parameterization. This involves using machine-learning-based image analysis to read the geometry of liver lobules and zonation patterns of steatosis from histopathological images, in addition to laboratory data to provide initial and boundary values for solutes and the ODE models as well as information on the transplantation process, e.g., cold ischemia time. The simulation yields spatial evolutions of the considered phases and solutes within the liver lobules over time based on the given data. Such an in-silico model, which connects structure, perfusion, and function of the liver via the interaction between mechanical properties of the graft, hepatic perfusion, and the affected metabolism could facilitate clinical decision-making cf. [2].

[1] B. Christ et al. Frontiers in Physiology 8:00906 (2017)

[2] B. Christ et al. Frontiers in Physiology 12:733868 (2021)

[3] W. Ehlers. „Foundations of multiphasic and porous materials” in Porous Media: Theory, Experiments and Numerical Applications (2002).

[4] T. Ricken et al. BMMB 14:515-536 (2015)

Various research works from the biomechanics community focus on the inflation of pressurized tubular materials as this problem can be directly linked to the modeling of arterial tissue (see, e.g. Merodio & Ogden). In the modeling, the material is often taken to be hyperelastic, and the relation between loading deformation is given in terms of a strain energy density function. This function is often decomposed to account for the contribution of the different types of constitutes. The pressurization of the tubular structure is accompanied by various types of responses in both the stable and unstable parts of the inflation (see, e.g. Topol et al. 2022).

Biological soft tissue consists of cellular and non-cellular components. The relationship between stress and strain in the material evolves, and is determined by viscoelastic and tissue remodeling processes (see, e.g. Topol et al. (2021), Wineman & Pence (2021)). Due to the interplay of countless biological and chemical processes, the mechanical properties of tissue will be time-dependent. These processes may be initiated, stimulated, and mediated by different physical processes, for example in the form of mechanical stimuli in static or dynamic form (Stoffel et al., 2017).

This work on the inflation behavior of pressurized tubes under a time-dependent behavior of the material. It can be shown that the inflation of the material shows a large variety of responses that depend on various factors, that include the geometry of the considered problem (H. Topol et al. 2019)). A certain focus is given to the influence of the fiber natural configuration of the inflation behavior, which may differ from that of the embedding ground substance.

References
J. Merodio and R. W. Ogden. Extension, inflation and torsion of a residually stressed circular cylindrical tube. Continuum Mech. Thermodyn. 28:157–174 (2016)

M. Stoffel, W. Willenber, M. Azarnoosh, N. Fuhrmann-Nelles, B. Zhou, B. Markert. Towards bioreactor development with physiological motion control and its applications. Med. Eng. Phys 39: 106–112 (2017)

H. Topol, H. Demirkoparan, T. J. Pence. Morphoelastic fiber remodeling in pressurized thick-walled cylinders with application to soft tissue collagenous tubes. Eur. J. Mech. A/Solids 77: 103800 (2019)

H. Topol, H. Demirkopararn, T. J. Pence. Fibrillar Collagen: A Review of the Mechanical Modeling of Strain Mediated Enzymatic Turnover. Appl. Mech. Rev. 73: 050802, 2021

H. Topol, N. K. Jha, H. Demirkoparan, M. Stoffel, and J. Merodio. Bulging of inflated membranes made of fiber reinforced materials with different natural configurations. Eur. J. Mech. A/Solids 94: 104670 (2022)

A. Wineman and T. J. Pence. Fiber-reinforced composites: nonlinear elasticity and beyond. J. Eng. Math. 127: 30 (2021)

Chiral molecules form a plethora of liquid-crystalline phases, even in equilibrium. Liquid-crystalline phases formed by chiral and active agents have features that are even more counterintuitive. For instance, while in passive matter, chirality tends to cloak itself from the long-wavelength elastic and hydrodynamic properties, activity reveals its effect in cholesteric and chiral columnar phases, leading to unique forms of odd elastic behaviours. Specifically, cholesteric phases display a unique odd elastic force density tangential to the contours of the constant mean curvature of layer undulation. This ultimately leads to the formation of an antiferromagnetically organised, columnar, vortex-lattice array in extensile active cholesteric fluids. Columnar materials which are both polar and chiral surprisingly display two-dimensional odd elasticity, even though the system is three-dimensional. The interplay of this odd elasticity with three-dimensional Stokesian hydrodynamics leads to an oscillatory optical mode.
The interplay of chirality and activity also leads to liquid-crystalline phases without a passive analogue, such as time liquid crystals. In this phase, planar-chiral active elements — spinners — that tend to align with each other trade in rotation symmetry breaking for time-translation symmetry breaking. That is, in these spontaneously rotating aligned states, rotation symmetry is restored in a time-averaged sense. Therefore, while such states are aligned, they escape the celebrated Simha-Ramaswamy instability that plagues uniaxial active suspensions by spinning out of unstable configurations. In my talk, I will discuss both new phases of chiral active liquid crystals and their properties and the effect of activity on classical chiral liquid crystalline phases. I will specifically highlight how the interplay of chirality and activity allows the former to affect dynamic and static properties of ordered states to a much greater degree than in passive soft materials, as summarised above.
Given that most biomaterials are chiral, biological systems ranging in scales from subcellular to all the way to organisms are the natural domain for detecting chiral active liquid crystalline organisations and their effect. Indeed, some of the phases that I will describe in my talk were subsequently observed in in vitro experiments with properties consistent with my predictions. Furthermore, given chiral active liquid crystals have material properties that can, in principle, make them useful for engineering applications, a further important avenue of research will involve synthesising them artificially.

Whether one considers swarming insects, flocking birds, or bacterial colonies, collective motion arises from the coordination of individuals and entails the adjustment of their respective velocities. In particular, in close confinement, such as those encountered by dense cell populations during development or regeneration, collective migration can only arise coordinately. Yet, how individuals unify their velocities is often not understood. Focusing on a finite number of cells in circular confinements, we identify waves of polymerizing actin that function as a pacemaker governing the speed of individual cells. We show that the onset of collective motion coincides with the synchronization of the wave nucleation frequencies across the population. Employing a simpler and more readily accessible mechanical model system of active spheres, we identify the synchronization of the individuals' internal oscillators as one of the essential requirements to reach the corresponding collective state. The mechanical 'toy' experiment illustrates that the global synchronous state is achieved by nearest neighbor coupling. We suggest by analogy that local coupling and the synchronization of actin waves are essential for emergent, self-organized motion of cell collectives.

In this talk, I will describe two biologically inspired systems that can be described using the same geometrical Hamiltonian formalism. The first is ATP synthase proteins which rotate in a biological membrane. The second is swimming micro-organisms such as bacteria or algae confined to a 2D film. I will show that in both cases, the active systems self-assemble into distinct structural states - the rotating proteins rearrange into a hexagonal lattice, whereas the micro-swimmers evolve into sharp lines with a particular tilt. While the two systems differ both on the microscopic, local interaction, as well as the emerging, global structure, I will show that their dynamics originate from similar geometrical conservation laws dictated by a Hamiltonian formalism applicable to a broad class of fluid flows.

Many active biological particles, such as swimming microorganisms or motor-proteins, do work on their environment by going though a periodic sequence of shapes. Interactions between particles can lead to the phase-synchronization of their duty cycles. We consider collective dynamics in a suspension of
such active particles coupled through hydrodynamics. We demonstrate that the emergent non-equilibrium states feature stationary patterned flows and robust unidirectional pumping states under confinement. Moreover the phase-synchronized state of the suspension exhibits spatially robust chimera patterns in which synchronized and phase-isotropic regions coexist within the same system. These findings demonstrate a new route to pattern formation and could guide the design of new active materials. An extension of the same theory for treating ciliated surfaces quantitatively captures the instabilities and flow pupmping behaviour of ciliated carpets and metachronal waves.

Axolotl salamanders regrow entire limbs throughout life using a molecular machinery similar to that used in the embryonic development of our limbs. A series of genetically-regulated molecular markers, known as morphogens, coordinate key mechanisms of this process. But exactly how morphogen patterning drives skeletal joint formation is not fully understood yet. Computational models provide a means to explore the factors influencing this process. In addition, we have access to three-dimensional (3D) microscopy images of morphogen expression in regenerating axolotl forelimbs obtained using a novel technique [1], which can help inform such models. But to effectively abstract key concepts and test hypotheses about how biological processes work using computational models, we must have a thorough grasp on how model parameters and conditions influence predicted outcomes.

One of the prevailing theories to explain the formation of morphogen patterns was proposed by the mathematician Alan Turing [2]. Starting from a nearly uniform initial state, and through suitable nonlinear interactions between reacting and diffusing morphogens, stable spatial patterns, known as Turing patterns, are obtained. Turing patterns have been used to investigate the formation of skeletal limb structures, but previous studies only used one or two-dimensional models, which have limited applicability due to the joint’s asymmetrical structure. Extending to 3D is challenging due to the wider variety of patterns and increased complexity.

We use linear stability analysis and finite element modeling to predict Turing pattern emergence in a 3D generic domain using the Schnakenberg activator-substrate model [3]. We provide a framework to identify the critical factors necessary for specific morphogen pattern emergence, and explore the role of initial conditions, model parameters and domain size on the predictions. We have observed that initial conditions on the activator have a stronger impact on the final pattern than initial conditions on the substrate, which has important implications for modeling morphogen expression in joint formation.

Our findings establish the groundwork for upcoming research, where we aim to employ the model developed here to predict the morphogen expression resulting from the axolotl joint formation experiments. Only through a comprehensive understanding of the factors influencing pattern emergence will we be able to successfully use Turing patterns as a model for morphogen patterning in joint formation.

References

[1] Lovely AM et al. HCR-FISH in Ambystoma mexicanum Tissue. In: Salamanders Methods in Molecular Biology. New York, 2023. p.109–22, doi:10.1007/978-1-0716-2659-7_6.

[2] Turing A. The chemical basis of morphogenesis. Philos Trans R Soc London, 1952; 237(641):37-72, doi:10.1007/BF02459572.

[3] Ben Tahar S at al. Turing pattern prediction in three-dimensional domains: the role of initial conditions and growth. Preprint on bioRxiv, 2023; doi: 10.1101/2023.03.29.534782.

Cells are constantly interacting with their environment, adapting their behavior in response to different stimuli and environmental conditions. This cellular adaptation occurs via changes in gene expression derived from changes in the physiological environment, giving rise to phenotypic plasticity. Phenotypic plasticity plays a key role in many steps of cancer progression, to the point that it has been included among the Hallmarks of Cancer. Indeed, it partly explains some of the most characteristic features of cancer, such as metastasis or drug resistance, probably the main challenges for improving cancer prognosis. Hence, understanding the mechanisms that trigger cellular adaptation is crucial for advancing in the study of this disease and its eventual treatments

In this line, mathematical models and simulation have great potential for gaining insight into complex cell process such as adaptation, and testing hypotheses regarding the effects of different environmental conditions in the adaptive response of tumors. In this work, we present a novel mathematical framework to model cellular adaptation and phenotypic heterogeneity in cell populations interacting with their environment. The model is based on the concept of internal variables, which are used to model cell state and regulate cell behaviour, in addition to Partial Differential Equations (PDEs) to describe the evolution of cell populations and the spatial concentration of chemical species from a continuum point of view. The proposed approach allows to consider not only cell response to environmental changes, but also reversibility and inheritance typical of phenotypic changes.

After presentation of the model, we particularize it to the case of glioblastoma (GBM) adaptation to hypoxia. GBM is the most common and lethal primary brain tumor. It has a dismal prognosis with a 5-year survival rate of 5%. The poor response of GBM to treatment is a consequence of its intrinsic and acquired drug resistance. This resistance is enhanced by hypoxia, a defining feature of GBM. Therefore, it is important to study how hypoxia governs cellular adaptation in GBM to improve our understanding on treatment response and, eventually, GBM prognosis.

The objective of this study is analyzing the potential of the derived model for capturing important biological trends in GBM evolution and adaptation under variable oxygen concentrations. After an extensive parametric analysis, different oxygenation conditions are tested. The results show the flexibility of the model for capturing the variability existing among GBM tumors. The model is also able to capture some observed experimental trends, such as the increased aggressiveness and resilience of tumors undergoing cyclic hypoxia.

In short, the developed framework presents an alternative to modelling cellular adaptation and may, with suitable validation, help in designing predictive tools as well as in silico clinical trials.

Last few years have brought into broader scientific audience examples of very robust and prominent periodic spatio-temporal patterns in biological systems [1,2]. Until recently, patterns such as spirals and planar wave trains were observed almost exclusively in pure chemical nonlinear systems. Reaction-diffusion equations are commonly used for understanding such phenomena. In these systems, onset of periodic wave trains is described by the so called wave instability (WI). Here we present generic three component activator-depleted substrate-inhibitor model which is a minimal model for patterns created by small GTPase signalling molecule RhoA. The model reproduces patterns which are observed experimentally [1] and explains some quantitative relationships between amplitude, temporal period and concentrations of some components.

Activity of RhoA may also induce contraction of cell surface (cortex) by activation of a motor protein myosin, which may result in advection of chemical components. Hence, our system becomes a reaction-diffusion-advection system with temporally and spatially dependent velocity field. We are interested in modification of chemical instabilities by biomechanical modes and interactions of those modes in simple reaction-diffusion-active gel framework. In this approach cell cortex is treated as a visco-elastic fluid in which active stress is controlled locally by the chemical subsystem [3]. We found that biomechanical coupling may increase parameter range in which spatially homogenous solution is unstable and periodic spatial patterns are possible. We show that interactions between chemical and biomechanical modes may lead to the change of spatiotemporal characteristics of patterns created by the chemical subsystem, their destabilization as well as formation of new patterns. Moreover, we discuss the role of the delay between RhoA activity and active stress. 2D simulations on the plane and stability analysis of 1D homogenous and nonhomogeneous solutions are presented.

[1] A. Michaud, M. Leda, Z. Swider, S. Kim, J. He, J. Valley, J. Huisken, J. Landino, A. Goryachev, G. von Dassow, W. Bement, A versatile cortical pattern-forming circuit based on Rho, F-actin, Ect2, and RGA-3/4. J. Cell Biol., 221(8):e202203017, (2022)

[2] W. Bement, M. Leda, A. Moe, A. Kita, M. Larson, A. Golding, C. Pfeuti, K-C. Su, A. Miller, A. Goryachev, G. von Dassow, Activator-inhibitor coupling between Rho signalling and actin assembly makes the cell cortex an excitable medium, Nature Cell Biol., 17(11), 1471 – 1483 (2015).

[3] Frank Jülicher et al, Hydrodynamic theory of active matter, Rep. Prog. Phys. 81, 076601 (2018).

Red blood cells (RBCs) are known as important formed elements in the blood. The functions of RBCs are to transport oxygen from the lungs to the tissues and metabolic waste, carbon dioxide, from the tissues to the lungs and to maintain systemic acid/base equilibrium. Along with that, RBC is known to release adenosine triphosphate (ATP) when it is subjected to shear stress. Subsequently, ATP molecules bind to purinergic receptors, to activate a cascade of biochemical reactions in endothelial cells (ECs) to release sequestrated ubiquitous calcium ion from endoplasmic reticulum (ER). The response of EC to ATP is in the form of transient calcium. Calcium is well known to regulate the activity of many enzymes in order to maintain the cellular homeostasis. Nevertheless, unregulated free calcium concentration could lead to serious pathological conditions such as cell death. In order to avoid such consequences, EC manages to maintain its physiological calcium concentration in the presence of ATP using ATPdriven pumps present in the plasma/ER membrane and the desensitization of the purinergic receptors. In order to understand how ATP released from RBCs affect the intracellular homeostasis in a vascular wall, we firstly developed a minimal calcium model, which guarantees the intracellular calcium homeostasis. Secondly, we couple it to RBC flow in a two-dimensional channel for a given imposed parabolic flow. In simulation, we use immersed boundary-lattice Boltzmann method (IB-LBM) to solve the fluid flow and the ATP release from RBC. We carried out several simulations varying flow strength, channel width, and concentrations of RBC (hematocrit) in order to emulate the blood flow in microcirculation. The endothelium helps maintaining a steady ATP concentration, avoiding the abnormal rise in the ATP concentration released from RBCs. With varying flow strength and hematocrit for a given channel width, we found that the ATP concentration and the cytoplasmic transient concentration increase with increase in the flow strength as well as hematocrit, and this leads to the cytoplasmic transient time decrease. Due to the relatively small peak times and amplitudes of cytoplasmic calcium at high flow strength as compared to that at low flow strength for all hematocrit, there is a possibility of calcium propagation from the high flow strength region to the low flow strength region for the coordination of cellular functions. Similarly, we observed a possibility of calcium propagation from low confined channels to the medium or highly confined channels for all hematocrit for a given flow strength. It would be interesting to carry out a further study in a vascular network in order to get more insights on the calcium propagation as each branch of vascular network may have unequal concentrations of RBC and the flow strength.

Spider silk is an extraordinary protein material that exhibits counterintuitive mechanical behaviors such as a reduction in stiffness of several orders of magnitude, supercontraction (i.e. a shortening of up to ~60% in length), and twist upon exposure to high humidity. These non-trivial responses originate from a unique polymeric structure made of crystalline domains that are embedded in a highly aligned amorphous matrix. Broadly, high humidity leads to water uptake by the silk, which motivates the dissociation of intermolecular hydrogen cross-linking bonds. In this talk, I will present energetically motivated models that explain the origin of supercontraction and twist in spider silk. Using tools from statistical mechanics, I will show that the dissociation of intermolecular bonds gives rise to a transition from a glassy to a rubbery phase, an increase in entropy, and a decrease in free energy. These factors shed light on the underlying mechanisms that govern supercontraction and agree with experimental findings. In addition, I will employ a continuum-based framework to show that the twist behavior originates from helical features that exist in a glassy spider silk fiber. The merit of these works is two-fold: (1) they account for the microstructural evolution of spider silk in response to water uptake and (2) they provide a method to characterize the microstructural evolution of hydrogen-bond dominated networks. The insights from the presented models pave the way to the design of novel biomimetic fibers with non-trivial properties.

1) N. Cohen, “The underlying mechanisms behind the hydration-induced and mechanical response of spider silk”, Journal of the Mechanics and Physics of Solids, 172:105141, 2023.

2) N. Cohen and C.D. Eisenbach, “Humidity-Driven Supercontraction and Twist in Spider Silk”, Physical Review Letters, 128:098101, 2022.

3) N. Cohen, M. Levin, and C.D. Eisenbach, “On the origin of supercontraction in spider silk”, Biomacromolecules, 22:993-1000, 2021.

The study of the mechanical properties of viscoelastic composites has been of great interest due to their unique characteristics. Following the methodology proposed in [1,2], we study the effective properties of nonlinear viscoelastic heterogeneous materials. With this goal, we employ the asymptotic homogenisation technique to decouple the equilibrium equation into a cell and a homogenised problem. The theory developed in this work is specialised to the case of a strain energy density of Saint-Venant type, with the second Piola-Kirchhoff stress tensor also featuring a viscous contribution. Within this setting, we frame the general theory in the case of infinitesimal displacements to make use of the correspondence principle which results from the employment of the Laplace transform. This choice is also advantageous to avoid the numerical complications arising in a finite theory. Furthermore, it permits obtaining the classical cell and homogenised problems in linear viscoelasticity as a special case. We frame our analysis by considering the case of uniaxially fibre-reinforced composites and, taking inspiration from [3], we write short formulae for the effective coefficients associated with the antiplane problem. Finally, after selecting different constitutive models for the terms associated with the memory of the constituents, our results evidence that the approximations of the semi-analytical method converge rapidly and comparisons with data available in the literature show a good agreement. Further developments of this work aim to generalise the model using the covariant formulation of continuum mechanics and to include a broader analysis of different microstructural geometric arrangements of nonlinear viscoelastic composites. A further scope is to frame the general theory in biological scenarios of interest. These include but are not limited to, biological fibrous tissues, such as muscles and connective tissue. It is expected that further research in this area will lead to new research questions in materials science and biomathematics.

[1] Pruchnicki, E. Hyperelastic homogenized law for reinforced elastomer at finite strain with edge effects. Acta Mechanica 1998, 129, 139–162. https://doi.org/10.1007/bf01176742.

[2] Ramírez-Torres, A.; Di Stefano, S.; Grillo, A.; Rodríguez-Ramos, R.; Merodio, J.; Penta, R. An asymptotic homogenization approach to the microstructural evolution of heterogeneous media. International Journal of Non-Linear Mechanics 2018, 106, 245–257. https://doi.org/10.1016/j.ijnonlinmec.2018.06.012

[3] Rodríguez-Ramos, R.; Otero, J.; Cruz-González, O.; Guinovart-Díaz, R.; Bravo-Castillero, J.; Sabina, F.; Padilla, P.; Lebon, F.; Sevostianov, I. Computation of the relaxation effective moduli for fibrous viscoelastic composites using the asymptotic homogenization method. International Journal of Solids and Structures 2020, 190, 281–290. https://doi.org/10.1016/j.ijsolstr.2019.11.014.

Stress shielding minimization is one of the major issues in implant design. Currently, aseptic loosening brought on by stress shielding is one of the primary causes of revision surgery. As bone regeneration is triggered by a stress stimulus, poor load transmission to the bone can result in a low stimulus, which can cause bone decay and other issues.

The development and evolution of additive manufacturing techniques have made it possible to significantly reduce stiffness by introducing porosity into implant materials as virtually any shape can be manufactured through those processes, regardless of shape complexity. In order to encourage cell adhesion and proliferation, porous geometries are frequently used in the construction of scaffolds. Additionally, porous geometries allow for the fine tuning of mechanical properties through changes to its topological design.

Feed-forward neural networks are able to do complex non-linear mapping between the input and output data. It has been shown that a neural network with one hidden layer and a number n of neurons in capable of representing any function. When a neural network presents several hidden layers, it is usually considered a deep learning approach.

The objective of this study is to create an optimal design by training a neural network to produce the ideal unit cell topology for a given constitutive elastic matrix. As a result, the network has the ability to reverse the homogenization process my mapping the relationship between the constitutive elastic properties and the unit cell geometry. A feed-forward neural network was created and trained in MATLAB with data generated from a set of several different geometries.

To each of these geometries, homogenization with periodic boundary conditions was performed. The lattice was modelled as a biphasic material where the solid phase was modeled with the properties of the material and the remainder area of the representative volume element (RVE) was considered to be a void phase with 1e-06 of the Young’s modulus of the solid material in order to reduce the influence of these elements to the homogenized constitutive matrix. A uniform mesh of square 2D elements allows to directly impose the periodic boundary conditions. The original geometry is therefore simplified to fit the uniform mesh where each element, a structured square, is either attributed solid or void properties.

The constitutive matrix used as the input of the network was obtained by applying a deformation gradient leading to a strain state where all components, but one is null, and thus, the macro-stress tensor of the RVE is equal to one line in the constitutive matrix. The linear-elastic analysis is run using ABAQUS as the solver.

Composite or microstructured materials have been long since considered as important means to engineer and optimize mechanical properties for specific applications. With the advent of additive manufacturing (also known as 3D-printing), production of artificial constructs conceived to possess specific optimal properties is now becoming possible, with applications ranging from construction to biomimetic materials. The architecture of such composites is typically based on designing features at a small scale, that lead to the desired large-scale behavior of structural elements. In light of this, theoretical studies of composites often involve asymptotic homogenization or alternative upscaling techniques. Most often, a periodic assembly of basic units is a practical way to build metamaterials.

In the design and analysis of composite materials based on periodic arrangements of sub-units it is of paramount importance to control the emergent material symmetry in relation to the elastic response. In numerous applications it would be useful to obtain effectively isotropic materials. While these typically emerge from a random microstructure, it is not obvious how to achieve isotropy with a periodic order. I has been long since recognized that inclusions distributed on a two-dimensional hexagonal lattice orthogonally extruded in the third dimension can give rise to transversely isotropic materials. Nevertheless, in spite of some computational evidence emerged in recent years, a generalization of this result to full three-dimensional isotropy has so far remained elusive.

We present a rigorous and yet simple proof of the fact that a periodic arrangement on a face-centered cubic lattice of spherical inclusions of an isotropic solid within an isotropic matrix gives rise to a large-scale isotropic response [Mech. Res. Commun. 126 (2022) 104001]. In doing so, we also show that any rhomboidal computational cell that generates such a lattice can be used to successfully design homogenized solids in which the material symmetry is not affected by the periodicity of the construction, since the latter would preserve even the largest possible symmetry group. It is significant to observe that the geometric symmetry group of such rhomboidal cells is strictly smaller than the symmetry group of the lattice they generate, but the lattice and not the cell is the geometrically relevant structure when analyzing large-scale properties.

We frame our discussion in the context of linear elasticity by introducing a normalized Voigt representation of the elasticity tensor which is very convenient for the identification of material parameters and symmetries for inclusion lattices. Such a representation hinges on the definition of a basis for the space of symmetric tensors that is linked to the generators of the periodic lattice. In this way, the coefficients of the 6 by 6 elasticity matrix that describes the linear stress-strain relation acquire a material meaning which is independent of the coordinate basis chosen to represent the tensors.

There is a shortage of physics-based tools, available to orthopaedic surgeons, to quantify their everyday decision-making measures towards resolving their patients’ orthopaedic disorders. Such decisions rely purely on static medical imaging and surgical experience gathered over the years. As such, there is general consensus about the lack of understanding surrounding the consequences of performed surgical interventions on the resulting patients’ biomechanics. Such problems are prevalent across orthopaedics, esp. with implants, where the influence of implantation on individual biomechanics is unknown. Hence, at this juncture, the subject-specific interaction of the implant with the patient’s locomotor system cannot be determined. Our motivation is to overcome these issues using high-fidelity, physiologically realistic, in-silico analysis of the patient’s musculoskeletal disorder at a given biomechanical joint, thereby enabling implant testing for an individual patient or even performing large scale clinical trials for a cohort of virtual patients. In this regard, I will be present our current work on Fraunhofer IPA’s In-silico Human Modelling platform (ISHM), ranging from medical imaging to complex 3d biomechanical simulations of the musculoskeletal system simulations with very detailed physiologically human models.

The biomechanics of every single joint in the human body is extremely complex. Their physiological motions within the framework of stable joint mechanics are balanced by a complex system of pre-stretched muscles, tendons, ligaments and various other connective tissues and controlled by targeted muscle activation. Such a biomechanical joint system reacts very sensitively to changes in the properties of its components. Therefore, forced joint models commonly used in simulations cannot represent the real joint biomechanics as they end up generating unphysical coercive forces. As a result, we get inaccurate muscle forces and joint kinematics, far from absolute physiological reality. In the current study, I will present some physiological joint FE analyzes of the foot, knee, and hip without joint constraints. The muscle-driven forward analysis, together with the physiological musculoskeletal model, allows us to understand joint mechanics much better, especially when it comes to the biomechanical analysis of subjects with musculoskeletal disorders. Patient-specific model generation, the biomechanics of soft- and hard-tissues and accurate representation of joint biomechanics are of significant importance for realistic FEA. Such analysis would pave the way for in silico engineering of medical products in orthopaedics such as implants, prosthetics and orthoses. Improved functional fitting of the products will enhance their performance, which shall have a positive impact on the patient’s biomechanics and their overall quality of life.

Bone tissue possesses a remarkable capacity for adaptation to external loads, allowing it to modify its structure and density accordingly. Cancellous bone, the spongy network of trabeculae that constitutes the inner part of bones, undergoes microstructural changes in response to under- or overloading, which may result in the reinforcement or narrowing of its constituent trabeculae and alteration of the microstructural pattern. Due to the microstructure, size-dependent effects may play a role. Multi-scale approaches consider bone as continuum matter and resolve the trabecular structure directly at the subscale. This is, however, algorithmically cumbersome to implement and, above all, involves additional computational effort.

Here we present a micromorphic approach that accounts for both the heterogeneous substructure of the material, without resolving it explicitly, and the nonlocality of the bone remodelling process, which is physiologically motivated by spatially correlated mechanosensing and regulation. Our approach has the advantage of being able to dispense with laborious neighborhood sampling, as is the case with integral approaches, and higher continuity requirements, as is the case with higher gradient approaches. Our approach is implemented in the open source finite element environment deal.II.

Our methodology employs nominal bone density as a macroscopic quantity for the bone mass to volume ratio in the underlying trabecular microstructure. This approach allows us to account for the heterogeneous nature of bone while avoiding the need to resolve individual trabeculae. In this model, we use the theory of open-system thermodynamics, which assumes a mass source proportional to the change in nominal density over time. This mass source is equated with a mechanical stimulus, and the stored energy is compared to an attractor, which represents a biological stimulus that drives remodelling. In the local case, the stored energy is a local quantity that depends on macroscopic deformation. In our novel non-local approach, we extend this to include a micromorphic and a scale-bridging component, allowing us to incorporate non-locality with a characteristic length scale for the heterogeneous microstructure and a scale-bridging parameter that couples the micro and macro deformation. This approach enables us to account for the interaction between continuum points and to model how points in the material that are not directly loaded react to the loading of their neighbors.

We present this approach in detail and demonstrate its efficacy using benchmark examples. We then apply this approach to long tubular bones and compare it with a series of CT images of femoral heads.

Mathematical models of the heart have evolved from single-physics representations on simplified geometries to coupled multi-physics representations of the whole heart with high anatomical fidelity. Here we present a fully-coupled electromechanical model of the human heart including all four chambers. Electrophysiology, active continuum biomechanics, and closed-loop circulation are modeled in a multi-scale approach. State-of-the-art models based on human physiology are used to describe membrane kinetics, excitation-contraction coupling and active tension generation in the atria and the ventricles. The validity of the model is demonstrated through simulations on a personalized whole heart geometry based on magnetic resonance imaging data of a healthy volunteer.

The proposed framework for the fully coupled cardiac electro-mechanical problem comprises a detailed description of appropriate boundary conditions such as a lumped parameter model of the human circulatory system and a contact handling that replicates the effects of the tissue surrounding the heart. To solve the coupled electro-mechanical problem, we apply a staggered scheme in time where the monodomain equation describing cardiac electrophysiology and the non-linear deformation resulting from the balance of active tension, passive material forces, chamber pressure and contact to the surrounding tissue are solved sequentially. Additionally, the proposed electro-mechanical whole-heart model is coupled to a 0D closed-loop model of the cardiovascular system. Here, we use a quasi Newton method to update the pressure values in all four chambers and reach convergence in fewer iterations compared to standard Newton methods.

We provide parameterizations for the fully-coupled excitation-contraction model for cells of the atrial and ventricular myocardium. Both the intracellular calcium transient and the tension development match data of human tissue preparations from literature. Coupling the 0D lumped parameter model of human circulation to all four chambers of the 3D electromechanical model enables a faithful reproduction of the major phases of the cardiac cycle as well as the characteristic figure of eight shape in the atrial pressure volume loops and flow patterns observed in clinical practice.
After introducing the coupled model, we provide application examples how such models can be used to generate mechanistic insight for clinical challenges. Selected examples include the in silico study of (i) the side effects of atrial ablation to treat electrical arrhythmia on whole-heart mechanical function and (ii) the electro-mechanical pathomechanisms underlying the breakdown of ventricular wringing rotation in heart failure. For (i), we provide biomechanical evidence that atrial ablation has acute effects not only on atrial contraction but also on ventricular performance. Therefore, the position and extent of ablation scars is not only important for the termination of arrhythmias but is also determining long-term pumping efficiency. For (ii), we show that isolated changes of the electromechanical activation sequence in the left ventricle are not sufficient to reproduce the rotation pattern changes observed in vivo and suggest that further patho-mechanisms are involved.

The crescent-shaped fibro-cartilaginous menisci located between the articulating surfaces of the femur and tibia are essential for the structural integrity of a healthy knee. Until recently, partial meniscectomy was considered the gold standard for treating meniscal lesions. However, due to the poor mid- and long-term outcomes of meniscus material dissection, surgical meniscus treatment paradigms have shifted in the last decade to promote healing through repair or regeneration. For cell colonization and meniscus replacement, various scaffolds have been proposed, including synthetic polymers, hydrogels, ECM components, and tissue-derived materials. So far, however, the outcomes of these studies seem inconsistent and indicate the need for a more fundamental understanding of the basic control mechanisms in cell-scaffold interactions under different environmental parameters.
In this context, we are working on simulations of a PDE-ODE system that model the dynamics of two cell populations involved in cartilage tissue production (adipose derived stem cells and chondrocytes). They are expected to migrate, proliferate, and differentiate within a perfusion chamber containing an artificial scaffold (poro-elastic medium) impregnated with a chemoattractant. To simulate this problem, we have used a first order discontinuous Galerkin scheme in space and an implicit Euler scheme in time.

This system is linked to a fluid problem, representing a mechanical stimulus. It is modeled by a Biot-Darcy system coupled to unsteady Stokes equations. To take into account all the coupled interface conditions, we have employed Nitsche's method.

In this macroscopic problem, several factors are important, most notably the stem cell differentiation that is expected to be mostly induced by mechanical stress inside the porous medium of the scaffold. Our goal is to identify all of the parameters of interest and investigate their impact. There are several methods available in this framework.
When dealing with discretized solutions of parameter-dependent PDEs, the sensitivities with respect to certain parameters of interest can be computed directly from the original problem, but each parameter of interest requires the solution of a new system: it is the "direct method". As a result, when dealing with a large number of parameters, the "adjoint method" may be a viable option.

Because the flow direction in our experiment alternates at a frequency of 1 Hz, the long-term simulation with a time span of up to 28 days is definitely a tremendous challenge. Thus, we will show how to efficiently combine reduced basis techniques and sensitivities computation with these methods to further reduce computational time. We will concentrate on non-intrusive reduced basis methods that do not require any changes to the High-Fidelity (HF) code. They only use the HF code as a "black box" solver. These adaptations will be numerically demonstrated with several results from our finite element model problem.

Embryonic development is a complex and fascinating process that has puzzled researchers for over a century. One of its greatest mysteries is how cells within a homogeneous mass can differentiate and organize into a wide variety of patterns. A key factor in this process is the role of morphogens. Morphogens are chemicals signals that cells use to communicate with each other, and they play a crucial role in determining cell fate and tissue specialization. Alan Turing proposed a model for how morphogens interact, known as the reaction-diffusion system [1]. This model relies on a diffusion-driven instability and nonlinear feedback between chemical species. Since Turing's time, researchers have used reaction-diffusion systems to model patterning in a variety of biological applications. In recent years, computational and experimental models have further validated the relevance and potential of this model.

The emergence of Turing patterns in a 3D growing domain has been sparsely investigated. The few studies exploring 3D Turing patterns indicate that the extension from 2D to 3D domains is not straightforward because the introduction of an additional dimension leads to a wider variety of patterns [2]. Previous studies do not focus on the biological context of embryogenesis and do not consider tissue growth. This gap in the literature highlights the need for further investigation on Turing patterns formation in 3D growing domains.

We explore the pattern evolution in 3D growing domains using finite elements analysis. As the growth process occurs on a much larger time scale than the reaction-diffusion one, we did not couple the domain growth with the pattern formation. Instead, for a given geometry, we computed the final steady-state pattern before growing the domain and relaunching the simulation. To model tissue growth two techniques were investigated: one simulating the apical tissue growth which is the cell formation at the distal or medial end and the second simulating anisotropic expansion. In the first case, new elements are added to one side of the mesh. In the second case, the whole mesh is stretched in one direction. Because these two different ways produce different pattern outcomes, pattern selection and emergence of bifurcations are affected due to the non-linearities of the system.

Overall, our study provides new insights into the role of different growth types when modelling 3D Turing patterns.

References:

[1] Turing AM. The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society of London. doi: 10.1007/BF02459572

[2] Ben Tahar S at al. Turing pattern prediction in three-dimensional domains: the role of initial conditions and growth. Preprint on bioRxiv, 2023; doi: 10.1101/2023.03.29.534782.

The consistent mathematical modeling of (active) biomechanical systems
at large strains based on partial differential equations is a challenging task.
In particular, the evolution of most biological systems is not just driven by mechanics, but is typically the result of coupling with other physical or chemical fields. Thus, the consistent coupling e.g. to reaction-diffusion equations or phase-field models becomes necessary to develop a general theory of out-of-equilibrium nonlinear thermodynamics for biological systems.

In this talk, we demonstrate that for dissipative biomaterials, coupled balance equations can be expressed as GENERIC or gradient systems. These frameworks describe the evolution of a system through thermodynamical potentials and geometric structures, encoding the reversible or irreversible nature processes using internal variables, such as concentration of biomolecules, phase-fields and growth variables, in addition to the elastic deformation. This point-of-view enables us to derive thermodynamically consistent coupled field equations systematically and reveals underlying coupling mechanisms. Moreover, the additional variational structure of the evolution equations can be exploited to establish structure-preserving numerical discretization schemes.

We discuss the applicability of the GENERIC framework and gradient systems for a biomechanical model describing the spatio-temporal evolution of brain atrophy in Alzheimer’s disease recently proposed by Schäfer, Weickenmeier, and Kuhl.
The model combines an anisotropic reaction–diffusion model for the time-dependent local concentration of misfolded tau protein with a large deformation shrinking model to predict the concentration-dependent regional tissue atrophy.
Similar models are relevant for the description of nonlinearly coupled transport and phase-transition processes in hydrogels and are therefore also highly relevant for biomedical applications.

The expression "Hierarchical fibre bundles" indicates structures with one principal direction. The structure is composed of multiple fibrillar elements, i.e. with one principal direction. Interfaces may or may not be present between the fibrillar elements. Finally, each fibrillar element is itself defined as a fibre bundle, made of fibrils with different mechanical and geometrical properties. The properties of each scale arise from the properties of the lower scale.
Hierarchical fibre bundles are abundant in nature across a variety of anisotropic biological materials, such as tendons, muscles and silk. These structures often exhibit exceptional fracture properties with outstanding toughness, thanks to the presence of structural proteins and/or to bio-mineralisation and their arrangement across different scales.
Analysing the fracture behaviour of complex fibre-based structures is still a challenging task, even when using advanced computational tools. Current numerical models often require a set of assumptions and (over-)simplifications. In particular, many models require the knowledge a priori of the fracture path. Biological materials, however, are characterised by complex internal structures with random features; moreover, for nano- and micro-bundles, it is often impossible to capture the fracture phenomenon empirically. The fracture path is thus often unknown, and it is one of the open questions which needs to be addressed using mechanical models. Here, we propose an alternative solution to model the fracture mechanics of fibre-based materials, which combines the mathematical Fibre Bundle Model (FBM) and Finite Element (FE) methods.
The Fibre Bundle model is a general theory of failure based on the statistical distribution of the critical elongations of the fibres forming the bundle. FBM assumes the bundle as a system of parallel one-dimensional elements, being clamped at their extremities and undergoing external forces along their main direction. Being a mathematical tool, the FBM is based on analytical formulations and thus it can be computed at practically instantaneous speed using computational solvers. In this work, we use the FBM to obtain general properties of the fracture behaviour of a fibre bundle undergoing elongation along its main axis, and to describe the stress-strain curves characterising the fibrils in the structure. The non-linear stress-strain curve obtained is utilised as an input for an FE approach, which aids us in comprehending the deformation and fracture behaviour of the fibrils and the interaction between different fibrils and the eventual surrounding interfaces.
We demonstrate the utility of the approach by applying it to spider-silk threads and bio-polymeric fibre bundles that are known for their non-linear properties. The enhanced numerical approach to the fracture mechanics of mixed fibre bundles will improve our understanding of the function of biological structures and bio-inspired materials across scales.

Constant demand for improvement of cardiovascular treatments and necessity for better understanding of soft tissue behaviour is frequently tackled with finite element (FE) modelling. The accuracy of FE simulations considerably depends on precise and reliable material parameters. They are obtained through the process of material characterization that consists of experimental testing and subsequent data analysis. When it comes to soft tissue, planar biaxial tensile tests are widely applied because they resemble in vivo conditions fairly well and material is assumed as anisotropic. In spite of planar biaxial tensile tests being the go-to type of experiment, they are still not standardized. The possible variations are numerous, from specimen size and different gripping mechanisms to the choice of preconditioning stretches, experimental protocols, and optimization procedures during data analysis. Limited availability of soft tissue requires simple and tiny square specimens. Before the start of experiment, specimens are mounted with hooks, rakes, or sutures and slightly preloaded to avoid the specimen hanging and to ensure that its orientation is perpendicular to loading axes. Application of the initial loads reflects in a different stress state compared to the reference configuration. This pre-stretching of the specimen cannot be measured due to the nature of experiment. Thus, measured initial force results in zero stretches during the material parameter estimation. However, these occurrences are mostly neglected. Either existing pre-stretches are disregarded or initial forces are annulled manually. Obviously, neither of the cases corresponds to real stress state.

In this work we developed the procedure to calculate mounting pre-stretches and consequently material parameters are determined more accurately. Standard optimization procedure is based on minimizing the least squares difference between experimental and model Cauchy stresses. Errors induced while calculating experimental Cauchy stresses lead to inaccurate material parameters. To avoid such errors, additional deformation gradient that represents mounted specimen in a state prior to experiment, has to be introduced. Since mounting pre-stretches are not measured, a new data fitting process is established with supplementary nested loop that will iteratively correct pre-stretches to satisfy equilibrium between initial experimental and model stresses. The procedures were first applied on the series of simulated virtual planar biaxial experiments where the exact material parameters can be set and compared to the obtained ones. Holzapfel-Gasser-Ogden strain energy function was used in virtual experiments while material parameters were based on the literature. Furthermore, data fitting procedures were applied to data gathered from biaxial experiments on aortic tissue samples. The analysis has shown significant difference between obtained material parameters, especially in cases where HGO structural parameters such as fibre family angle and in-plane dispersion are considered as fitting parameters. The rate of error increases with the amount of applied pre-stretches and also decreases with increase of maximum achieved experimental stretches. Finally, annulling initial force proved to be better approach if only standard optimization is to be applied.

The dura mater is the outermost layer of the meninges, a layered membrane that covers and protects the brain and spinal cord. It is a dense and fibrous connective tissue that is composed of collagen and elastin fibers, proteoglycans, and other extracellular matrix components [1]. The dura mater plays an important role in maintaining the stability and integrity of the brain and spinal cord, as well as in regulating the cerebrospinal fluid dynamics and protecting the neural tissue from mechanical stress and injury [2]. Despite its important role in, e.g., traumatic brain injury pathology it is frequently neglected in computational and physical human head models. Until now, the biomechanical failure tests have not yielded information about the anisotropic mechanical behavior of various locations of the human dura mater under physiological conditions. Additionally, limited data exists on the orientation of collagen fibers, which has not been thoroughly quantified.

Therefore, this study aimed to investigate the mechanical and structural anisotropy of different locations of the human dura mater and to correlate mechanical data with the quantification of the orientation of collagen fibers. To achieve this objective, sixty samples from six donors were subjected to quasi-static, uniaxial extension tests until failure in a heated tissue bath, utilizing a Z020 torsion multi-axis testing system (ZwickRoell AG, Ulm, Germany) together with an Aramis image correlation system (GOM, Braunschweig, Germany). Additionally, the collagen microstructure of samples from four donors was analyzed using second-harmonic generation imaging.

The data obtained were used to determine the failure stress and strain, E-moduli, and a microstructurally motivated material model was employed to examine local differences in both structure and mechanics. Among others, significant differences in both collagen fiber dispersion and main fiber orientation were found. The structural parameters obtained from only four out of six donors yielded good fitting results to the structurally motivated mechanical model for the remaining two donors in all directions. This study establishes a foundation for further research on the microstructure and mechanical properties of the human dura mater, enabling realistic modeling and prediction of tissue failure.

References

[1] J. T. Maikos, R. A. Elias, and D. I. Shreiber, "Mechanical properties of dura mater from the rat brain and spinal cord," J Neurotrauma, vol. 25, no. 1, pp. 38-51, Jan 2008, doi: 10.1089/neu.2007.0348.

[2] D. J. Patin, E. C. Eckstein, K. Harum, and V. S. Pallares, "Anatomic and biomechanical properties of human lumbar dura mater," Anesth Analg, vol. 76, no. 3, pp. 535-40, Mar 1993, doi: 10.1213/00000539-199303000-00014

Cardiovascular diseases are still one of the major cause of death worldwide. They are responsible for significant morbidity and therefore both noninvasive and especially invasive treatments are still being developed. Our research group is working in the field of new materials for cardiovascular surgery.

There are several procedures in the field of cardiovascular surgery which are based on the mechanical interaction of the arterial wall and the external support. A typical example is Pulmonary Artery Banding (PAB) in infant patients. PAB is a palliative procedure that reduces pulmonary over-circulation in neonates suffering from certain congenital heart diseases and constitutes the first stage of intervention prior to the complete repair of cardiac defects. The principle of PAB consists in the implantation of a band around the pulmonary artery, which results in the reduction of blood flow.

Our goal is to develop resorbable fabric based on composite material consisting of a polylactide and polycaprolactone (PLA/PCL) copolymer nanofibrous reinforcement combined with a collagen matrix for use in cardiovascular surgery as the pulmonary artery banding material. The resulting composite composed of these polymers prepared by electrospinning followed by the impregnation of the layers with collagen (type I) dispersion is then collagen-poly[L-lactide/caprolactone] (COLL-PLCL). This compound is bioresorbable and gives the advantage of PAB as a one-step procedure.

In our work, we present the results of the comparison of the mechanical properties of our COLL-PLCL material with two kinds of commercially available Gore material that can be used in pulmonary artery banding. They are the Gore Pericardial Membrane and the Gore Cardiovascular Patch, and both materials are chemically based on the form of polytetrafluoroethylene (PTFE).

The strips of all three materials underwent uniaxial tensile tests. In the case of the COLL-PLCL material, it was a monotonic tensile test, whereas the Gore materials were tested cyclically. The COLL-PLCL materials were tested in the hydrated state.

The results of our experiments show that all three materials exhibit an approximately linear mechanical response under primary loading. Cyclic loading led to so called preconditioned state, in which the Gore materials exhibited some non-linearities in their response. The preconditioned Gore-Tex materials were also considerably stiffer than those in the pristine state. Gore materials showed higher stress at failure than COLL-PLCL material whereas strain at failure was approximately the same. What is significant is that the COLL-PLCL material was much more compliant than Gore. We hypothesize that it will therefore cause a less mechanobiological response to compliance mismatch during banding.

Although the biomechanics of the circulatory system is a very dynamically developing field, the study of venous aneurysms is somewhat outside the mainstream. Although aneurysm, a bulge of the blood vessel wall, is mainly known as a disease of the arteries, it can also be encountered in veins. It occurs mainly at the site of an arterio-venous shunt, which is artificially created as a so-called vascular access for hemodialysis. Other venous aneurysms are quite rare. However, the number of patients for dialysis is increasing, so the need to study tissue biomechanics in vascular access is a current scientific problem.

The exact mechanism of arterio-venous shunt aneurysm formation, apart from the rather vague claim that it is a consequence of altered hemodynamic conditions, is still unknown. The same applies to the constitutive models and their parameters found for such affected vein walls. However, it should be emphasized that not only the change in hemodynamics due to arterialization of the vein, but also the repeated needling of the tissue during the hemodialysis procedure are involved in the formation of the arteriovenous shunt aneurysm. These lead to significant scarring of the vessel wall, which results in a change in mechanical properties compared to healthy vein. As with any scar, we should expect a significant increase in connective tissue content.

The long-term goal of our study is to experimentally identify material models for pathologically affected veins and compare them with the behavior of healthy tissue. We then intend to use the constitutive models for the arterio-venous shunt aneurysm tissue in computational fluid dynamics simulations with a deformable wall to determine the changes that occur in the hemodynamics of the aneurysm at the vascular access.

As a first step in our study, we perform tensile experiments with vein wall samples, both healthy ones and veins with a created arteriovenous shunt. We chose to model the tissue behavior as elastic after preconditioning. The constitutive model used for the vein wall is adopted from the work of Gasser Holzapfel and Ogden (2006), or we model the vein wall as an incompressible, anisotropic, hyperelastic material that exhibits significant material nonlinearity. Because aneurysmal vein wall samples are collected during surgery from living patients, it is not possible to compare them with other donor veins, as they are not naturally available. Comparisons will be made against the characteristics of veins taken to create an aorto-coronary bypass graft. The number of experimentally measured samples in both groups is not yet large enough to achieve significant population differences. Therefore, our present results should be considered as preliminary.

Tissue fibrosis, classically associated with excessive accumulation of extracellular matrix (ECM) by activated fibroblasts, can occur in several pathological conditions. Fibroblasts are mechanosensitive cells and increased stiffness of the ECM activates fibrotic pathways [1]. To mimic the three-dimensional cell microenvironment in vitro, fibroblasts are cultured in 3D spheroids. Mechanical assessment of cell spheroids is accomplished via parallel plate compression and mechanical properties in a continuum sense are determined from force vs displacement (F-δ) curves. Current analysis methods are often limited to Hertzian theory and its modifications, that only account for small deformations [2]. However, spheroids exhibit low stiffness, thus undergoing large deformations at small external forces, and therefore current contact mechanics models fail to describe such behavior [3,4]. Here, we use the Tatara model, which considers both nonlinear elasticity and large deformations to describe the mechanics of cell spheroids under parallel-plate compression.

Cell spheroids consisting of primary human fibroblasts cultured for four days were subjected to parallel-plate compression testing (MicroSquisher, CellScale) fitted with a round tungsten cantilever and accompanying SquisherJoy V5.23 software (CellScale, Ontario, Canada). The fluid bath test chamber was filled with sterile phosphate buffered saline (pH=7.4). Stage and optics were calibrated according to manufacturer’s instructions. Samples were compressed up to 50% apparent linear strain at different displacement rates. F–δ data was fitted using linear least squares regression on the Hertz and Tatara model and its extended version (custom MATLAB code) with fully constrained contact points (F =0, δ =0). Images were captured via two digital cameras during compression and static states to determine contact radius and lateral expansion of cell spheroids over a range of deformations.

The dependence of compression modulus on the displacement rate was examined. The Hertzian (and extended) theory was successfully applied to the small strain regime of F-δ curves to extract the compression modulus. In contrast to the Hertzian theory, for compressive displacements over 10-25% (depending on displacement rate), where the force follows the third and fifth power of the displacement, the Tatara numerical analysis theory better described the non-linear behavior of cell spheroids. The predicted contact radius of cell spheroids was found to be in good agreement with the data obtained from image analysis. However, the lateral expansion was underestimated by the Tatara model, a difference which may arise from the adopted constant-volume assumption.

The Hertzian theory and its modifications can be applied for strain in the range of 10%-25% while, for larger deformations (above 25% strain), the extended Tatara model better describes the deformation behavior of cell spheroids also providing their lateral expansion. While Hertzian contact theory and its modifications is applicable for small deformations, hyperelastic and nonlinear elasticity theory should be employed to better describe the data at larger deformations.

1. M. Jones et al, eLife, 7 :1-24, 2018

2. Y. Efremov et al, Scientific Reports, 7:1-14, 2017.

3. Y Tatara et al, JSME INT A-SOLID M, 36:190-196, 1993.

4. K. Liu et al, J. Phys. D, 31 :294-303, 1998

Background. Despite biomechanical factors playing a pivotal role in the onset and development of disease [1], little is known concerning the fracture behaviour of most biological tissues. Most experimental testing techniques fail to acquire data for the mechanical fracture characterisation of soft tissues. Therefore, we propose the novel symmetry-constraint Compact Tension (symconCT) test for aortic tissue fracture testing. It allows for stable crack propagation and the in-depth exploration of fracture processes. Finite Element Method (FEM) modelling was combined with Digital Image Correlation (DIC) to identify the specimen-specific constitutive properties.

Material and Methods. The classical CT test was augmented with an elastic metal beam connecting the two clamps [2]. It pre-strained the 30x35 mm² specimens and forced the fracture to develop in the middle of the specimen. Specimens with a 10 mm pre-notch were obtained from the intima-media compound of the porcine abdominal aorta. Tensile loading, either in the axial or circumferential direction, was applied orthogonal to the pre-notch at 3 mm min^-1 (ADMET eXpert 4000 Universal Testing System). 2D FEM models were developed (Abaqus 2019) for six specimens (representing mean ± standard deviation of recorded symconCT data sets), incorporating a viscoelastic Yeoh model and an isotropic bilinear cohesive zone model. Whilst the viscous properties were fixed, the two elastic Yeoh coefficients, c₁ and c₂, were calibrated against the clamp force and DIC-based Green-Lagrange strains. The cohesive strength and the fracture energy were also calibrated against the force.

Results. Given the load aligned with the main collagen fibre direction in the media, circumferentially loaded specimens reached higher forces, with a mean of 5.76 ± 0.69 N against 3.89 ± 0.57 N of the axially loaded specimens (t-test, p<0.05). Whilst circumferential loading resulted in a zig-zag crack, axial loading provoked a straight fracture. The DIC showed the maximum principal strain being orthogonal to the crack tip. The identified elastic properties were c₁ = 34.70 ± 1.52 kPa and c₂ = 98.70 ± 5.63 kPa under circumferential loading and c₁ = 14.54 ± 5.24 kPa and c₂ = 67.34 ± 6.65 kPa under axial loading. The cohesive strength resulted in 235.62 ± 24.04 kPa and 242.40±72.48 kPa, the fracture energy in 1.57 ± 0.82 kJ m² and 0.96 ± 0.34 kJ m² under circumferential and axial loading, respectively. Given the calibrated FEM models, the FEM crack extensions closely matched the experimental ones.

Conclusions. The symconCT test allowed for stable crack propagation, permitting the exploration of the fracture behaviour ahead of the crack tip. The DIC-based principal strains revealed mode-I failure dominating the fracture in both loading conditions. The FEM models accurately simulated the tissue fracture. Although inertial effects were overaccelerated (mass-scaling), good agreement with the experimental recordings was observed. Additionally, the vessel wall shows anisotropic mechanical properties [3], characteristics which could have further improved our results.

References.

[1] N. Baeyens, M. A. Schwartz, Mol. Biol. Cell, vol. 27, pp. 7-11, 2016.

[2] M. Alloisio, M. Chatziefraimidou, J. Roy, T. Gasser, Acta Biomaterialia (under review).

[3] T.C. Gasser, Vascular Biomechanics, ISBN: 978-3-030-70965-5. Springer 2021.

The ability of bone to sense and react to its mechanical environment is well-known and undisputed. In qualitative terms, the prolonged exposure of bone tissue to increased mechanical loading (with respect to a "normal" physiological load level) leads to a corresponding increase in bone mass; usually associated with the filling of the bone tissue-surrounding pore spaces by additional bone tissue. A decreased mechanical loading, on the other hand, causes a corresponding depletion of bone tissue, leading to enlarged pore spaces and a decreased bone mass.

Several stimuli have been proposed as potentially relevant for the mechanobiological regulation of bone tissue, including direct cell stretching, hydrostatic pressure, and fluid flow-induced shear stresses. Furthermore, it is well known that bone tissue exhibits piezoelectric properties, and it has been suggested that piezoelectric excitation may be a contributing factor to the mechanobiological regulation of the activities of cells residing in the bone pore spaces as well. In this contribution, the focus is on the latter proposition. Due to the very small length scales at which the related processes occur, direct experimental validation or falsification of this hypothesis has turned out to be out of reach. Here, we investigate this potentially important mechanism by means of continuum micromechanics-inspired multiscale modeling.

In particular, the concept of continuum micromechanics, originally developed for elasticity upscaling (or, homogenization), was utilized for deriving and homogenizing a so-called electromechanical tensor, comprising stiffness as well as electrical quantities, and for up- and down-scaling mechanical and piezoelectric properties. This modeling concept was applied to bone tissue, spanning thereby the hierarchical organization from the molecular level of bone (where the so-called elementary constituents of bone can be identified) to macroscopic bone tissue. Applying this model to the multiscale micromechanical representation of bone tissue has revealed that the electrical stimuli arriving at the level of bone cells, in response to physiological macroscopic loading, is much too small for effectively stimulating bone cell activities. This result suggests that piezoelectric effects may be a contributing, but not a major factor for bone cell excitation. This contribution is concluded by an interpretation of the obtained results in the context of other mechanical stimuli - i.e., hydrostatic pressure and fluid flow-induced shear stresses - previously investigated also by means of continuum micromechanics-inspired multiscale methods, and by comparing their respective importance for the bone remodeling process (as suggested based on the computational results), revealing that the pressurization of bone pore spaces may be of bigger relevance than the other stimuli.

For accurate quantification of mineralised tissues using polychromatic lab x-ray sources, one needs to do two things: 1) to integrate for a very long time to get a good signal to noise ratio and 2) to calibrate the x-ray source and the detector response. At Queen Mary we build our own X-ray Micro-Tomography (XMT, synonym µCT) systems specifically designed for accurate quantification of mineral concentration. We build on the pioneering work of Professor James Elliott, the co-inventor of x-ray micro-tomography. We use CCD cameras and time delay integration to allow long integration times, and calibrate every scan with a multi-metal carousel, enabling us to correct beam hardening effects, model the x-ray emission spectrum and model our detector response which allow us to derive the sample Linear Attenuation Coefficient (LAC) and hence mineralisation.

Getting improved contrast [i.e., mineral concentration] resolution is crucial to all applications of XMT to studying skeletal tissues and especially where there are mixed tissue types with different degrees of mineralisation. The common clinically important fractures of the cortices of vertebral bodies and the femoral neck involve variable amounts of calcified ligament: this, and calcified cartilage in the end plates, reach higher levels of mineralisation than bone. They also have different histology, ‘grain’ and fracture properties. With our own systems, we can distinguish these tissue phases, but this is not possible in commercial systems. Mistakes will therefore be made is measuring cortical thickness – never mind that the tissue is not bone.

Context. Bone disorders like osteoporosis as well as active lifestyles increase the occurrence of bone fractures and joints damage. In orthopaedic and dental surgery, the implantation of biomaterials within bone tissue to restore the integrity of the treated organ has become a standard procedure. Their long-term stability relies on osseointegration phenomena, where bone grows onto and around metallic implants, creating a bone-implant interface. The bone tissue is a highly hierarchical material which evolves spatially and temporally during this healing phase. A deeper understanding of its biomechanical characteristics is needed, as they are determinant for the surgical success. In this context, we propose a multiscale homogenization model to compute the effective elastic properties of bone tissue (macroscale) as a function of the distance from the implant, based on the structure and composition at lower scales.

Methods. The model considers three scales: mineral foam, ultrastructure, and bone tissue. The elastic properties and the volume fraction of the elementary constituents of bone matrix (mineral, collagen and water), the orientation of the collagen fiber relatively to the implant surface and the mesoscale porosity constitute the input data of the model. At each scale, the continuum micromechanics theory based on the famous Eshelby’s representation of the uniform elastic field inside the ellipsoidal inclusion is applied. Experimental data were obtained from Ti6Al4V coin-shaped implants that were osseointegrated into cortical rabbit bones. The mineral platelet orientation -assumed to be parallel to the collagen fibers- and the mesoscale porosity were retrieved using small-angle X-ray scattering (SAXS) and using light microscopy, respectively (Le Cann et al., Acta Biomater, 116:391–399, 2020). The effect of their spatial variation on the bone anisotropic properties in the proximity of the implant were investigated.

Results and discussion. The findings revealed a strong variation of the components of the effective elasticity tensor of the bone as a function of the distance from the implant. The effective elasticity is primarily sensitive to the porosity (mesoscale) rather than to the collagen fibers orientation (submicroscale). However, the orientation of the fibers has a significant influence on the isotropy of the bone tissue, leading to a high degree of anisotropy when the fibers are oriented with an angle close to 45° with respect to the implant surface. When analyzing the symmetry properties of the effective elasticity tensor, the ratio between the isotropic and hexagonal components is determined by a combination of the porosity and the fibers orientation. A decrease in the porosity leads to a decrease of the bone isotropy and to an increase of the impact of the fibers’ orientation.

Conclusions. These results demonstrate that the collagen fiber orientation affects the effective elastic properties of the bone throughout the remodeling process in the proximity of an implant. Collagen fiber orientation should be taken into account to properly describe the effective elastic anisotropy of bone at the organ scale.

Magnesium implants appear as promising technology for load-carrying, bioresorbable bone regeneration tools (Hofstetter 2014 [1], Kraus 2012 [2]). In this context, quantitative exploration of the biomechanics and mechano-biology of the bone-implant compound system is of great interest. As a corresponding contribution, we here present an analytical beam theory representation (Pircher 2021 [3]) of a cylindrical magnesium implant located roughly in the middle of a rat femur and approximately orthogonal to the axis of the femoral shaft, which we map onto CT and SAXS-tomographic images of the pre-clinical animal test setting (Liebi 2021 [4], Grünewald [5] ). This reveals the correspondence between mechanical stress state distributions and micro-texture (collagen fiber orientation) across the bony organ. In more detail, the femur is first represented as a polygon with two straight lines, one associated with the femoral neck, and the other one associated with the femoral shaft. This polygon is the basis of a beam representation, where the shaft line is fixed at the knee, and the femoral head (the end of the neck) is loaded by external forces arising from standing. For the stress analysis, the implant is represented by the geometrical object "solid circular cylinder“, such that the cylinder’s volume, its radial normal component of the inertia tensor, its center of gravity, and its principal axes are identical to the corresponding quantities describing all the voxels showing the implant in the CT image. In the same way, the shaft region in proximity of the implant is represented by (i) a full cylinder representing the entire shaft, and (ii) a hollow cylinder representing the cortical shaft compartment only. Then, classical Bernoulli-Euler beam theory is applied to the bony portions of the implant-bone compound structure, with the beam cross-sections being dictated by the interaction of the roughly orthogonally positioned geometrical objects. As a measure for correspondence between stress states and texture directions, we evaluated the dot products between principcal stress directions and directions of maximum shear stress. As a result, somewhat surprisingly, not the principal stress directions, but the loading directions associated with maximum shear stress appear as the key mechano-biological drivers. Accordingly, the shear tractions of largest magnitude act orthogonal as well as parallel to the main texture direction, i.e. the collagen fiber direction.

[1] Hofstetter et al. (2014) JOM 66.4, pp. 566–572, 10.1007/s11837-014-0875-5.

[2] Kraus et al. (2012) Acta Biomaterialia 8.3, pp. 1230–1238, 10.1016/j.actbio.2011.11.008

[3] Pircher (2021) Thesis Wien, 20.500.12708/79654

[4] Liebi et al (2021) Acta Biomaterialia, 134, 804-817, 10.1016/j.actbio.2021.07.060

[5] Grünewald et al (2016) Biomaterials 76, pp. 250–260, 10.1016/j.biomaterials.2015.10.054.

Bone remodeling is a very complex and fine-tuned process, which is necessary to ensure a healthy bone structure. If this process gets out of balance – e.g., because of hormonal disbalance or the impact of bone metastases – pathologies like osteoporosis can appear. In this contribution we introduce a novel computational approach to investigate this balance by connecting the bone remodeling process with its microenvironment. Our goal is to better understand the well-balanced and complex dynamic of the subprocesses involved in healthy bone remodeling.

We implement a 3D stochastic cellular automaton (SCA), where voxels interact only with their nearest neighbors in a scaffold representing bone tissue. At each time point, each voxel can take one of four different states that stand for the different phases of bone remodeling: formation, quiescent bone, resorption, and environment. To create a compact representation of the frequency-dependent interaction of those voxel states we make use of methods borrowed from evolutionary game theory for the update rule of the cellular automaton [1]. This representation encodes knowledge about the mutual impact the different actors of bone remodeling (osteocytes, osteoclasts and osteoblasts) have on each other. Each parameter in the model has therefore a direct connection to the biological processes.

First, we set up simulations of the model with either only resorption or only formation. This choice reduced the model complexity and allowed us to determine parameter spaces for a self-regulating behavior for each of them. The self-regulating behavior is defined by resorption or formation starting and ending without further parameter tuning. Parameters outside the range of self-regulation will lead to either osteolytic lesions (resorption) or heterotopic ossification (formation). Further analyses supported the approach of a spatial model with a small neighborhood to simulate the local phenomena observed in bone remodeling.

Next, we coupled the two processes of resorption and formation. In the limit of separation of time scales, our model showed that self-regulating resorption followed by self-regulating formation reproduces the physiological bone remodeling behavior. Further analysis will create a more fluid coupling of the two processes while involving more parameters.

The model has the potential to use the role of the microenvironment to evaluate the impact of additional factors, such as drugs or bone metastases. We are planning on using experimental in vivo data from a breast cancer bone metastasis mouse model [2], which includes spatial and temporal dynamic of early osteolytic lesions, to fit additional parameters. Hopefully, these findings will add to the discussion, how pathological behavior might be controlled, if not even reversed.

[1] M. D. Ryser and K.A. Murgas, Bone remodeling as a spatial evolutionary game, Journal of Theoretical Biology, 2017

[2] S. A. E. Young, A.-D. Heller et al., From breast cancer cell homing to the onset of early bone metastasis: dynamic bone (re)modeling as a driver of osteolytic disease, bioRxiv preprint

The mechanical quality of trabecular bone is influenced by its inhomogeneous mineral content and spatial distribution at a microscopic scale. The bone remodelling process, which is the concerted action of osteoclastic bone resorption followed by osteoblastic bone formation, controls bone turnover. During the bone formation process (FP) the deposited organic collagenous matrix (i.e., osteoid) becomes mineralised. The latter process is regulated by the mineralisation kinetics which exhibits two distinct phases: a fast primary mineralisation phase lasting for several days to a few weeks and a secondary mineralisation phase that can last from several months to years. Variations in bone turnover and mineralisation kinetics can be observed in the bone mineral density distribution (BMDD), which can be used to distinguish healthy from pathological bone tissue at the scale of the bone matrix. Here, we propose a discrete statistical spatio-temporal bone remodelling model to study the effects of activation frequency (Ac.f) and mineralisation kinematics on the BMD distribution. In this model individual basic multicellular units (BMUs) are activated discretely on trabecular surfaces which then undergo typical bone remodelling sequence, i.e., resorption, reversal, osteoid formation, and mineralisation, with the latter process following a double exponential law. Our simulation results highlight that trabecular BMDD is strongly regulated by Ac.f and kinetics of secondary mineralization (t₂) in a coupled way. Indeed, Ca wt% increases with lower Ac.f and fast secondary mineralisation, while lower Ca wt% values are obtained for higher Ac.f and slower secondary mineralisation. Accordingly, the dynamic equilibrium can be achieved for different combinations of Ac.f and t₂. For example, a mean Ca wt% of 25, which has been reported in the literature based on qBEI experiments, can be obtained with Ac.f = 4 BMU/year/mm³ and t₂ = 8 years or with Ac.f = 6 BMU/year/mm³ and t₂ = 6 years. This close link between Ac.f and t₂ on BMDD results shows the importance of taking both characteristics into account in order to draw meaningful conclusions about bone quality. Indeed, pathological conditions such as osteoporosis demonstrate a similar pattern. We investigated post-menopausal and senile osteoporosis (type I and type II, respectively) and hypothesised that high Ac.f or very long formation period (FP) would result in type I osteoporosis, whereas underfilling and no filling would result in type II osteoporosis. Our results show that the apparent density and bone mineral fraction were similar for all osteoporosis hypotheses except for the non-filling one. However, when examining BMDD, significant variations were observed, especially in type I simulations. High Ac.f resulted in low Ca wt%, consistent with post-menopausal osteoporosis, while large FP led to high Ca wt%, which is only seen in type II osteoporosis.

The unique and somewhat enigmatic properties of biological apatites have made them an attractive area of research for many decades, with studies aiming to provide a comprehensive understanding of apatite physiochemistry and biological properties. Our work, however, initially began unconventionally at the end of the story, through investigating how the unique properties of biological apatites could be applied to the fields of forensic science, medicine and archaeology and finally ending with understanding the fundamental features of these properties. Our work with Professor John Clement began in the late 1990’s, with a collaboration between Professor Keith Rogers and Dr Sophie Beckett, on species differentiation of heated bone. The results of the study clearly demonstrated that apatite characteristics of unheated and heated bone exhibit significant inter-species variation, quantifiable using X-ray diffraction analysis. The work also highlighted the potential capability of distinguishing human from non-human bone based on apatite physicochemistry. Discussions with John throughout this project centred around ‘apatites not being a fixed piece of chemistry’ and there being ‘a need to understand the fundamental features of apatites to truly appreciate the applications of the work’.

Dr Charlene Greenwood’s PhD project followed, which aimed to provide a new insight into the fundamental mechanisms and processes associated with physicochemical changes to bone during heat treatment. The work, which studied both unheated and heated bone mineral chemistry, considered a controversial in vivo crystallite size control mechanism for biological apatites and mathematically (through Arrenhius and Johnson-Mehl-Avrami equations) described crystallisation kinetics of apatite during heating. Of importance from this work, was the role of carbonate in biological apatite formation and crystal growth, with Dr Emily Arnold’s PhD furthering our understanding of carbonates within the apatite lattice through the application of X-ray pair distribution functions. The work allowed the carbonate within the apatite lattice to be examined on a local scale during heating, with the results suggesting that carbonate does affect intermediate (and therefore likely average) order, but not local order. This insight is significant when considering apatite physiochemistry and disease, and area of research our team is particularly focused on.

Our work in the field of apatites and disease focuses on understanding apatite physicochemisty and disease initiation and progression. Our work on osteoporosis, which was heavily supported by John and the Melbourne Femur Research Collection, identified new fracture risk biomarkers based on tissue features, with differences in apatite crystal chemistry, nanostructure and microstructure identified between fracture and non-fracture groups. John’s legacy has continued with our work on ectopic calcification chemistry and its association with cancer. The team are currently exploring the view that calcification physicochemical characteristics could contain additional or independent source of diagnostic information, and we are applying this hypothesis to both breast and prostate microcalcifications to identify new chemical biomarkers for early diagnostics and prognostics.

Our presentation aims to discuss our fascinating journey with apatites and John, from the applications to their fundamental properties, with the work and John’s legacy continuing, even if we did take ‘the long way round’.

Introduction: Quantifying the effect of a new drug in the presence of natural variability between individuals is a major challenge of drug development. For bone diseases, where mechanics is an important outcome, one approach is to apply a non-intrusive uncertainty quantification method such as Monte Carlo on finite element (FE) models of bones, but with the drawback that a large number of simulations are usually necessary [1]. In this study, a stochastic FE (sFE) solver was developed that uses an intrusive method to predict the variation in bone mechanics given the variability in bone tissue elastic properties.

Method: The sFE solver is developed based on ParOSol, an open-source matrix-free FE solver specialized in bone mechanics. ParOSol solves the system of equations without storing the calculated matrix, resulting in a minimum memory footprint. Developing the sFE solver necessitates significant development in different modules from pre-processing to post-processing using C++ and Python programming languages. The core of the sFE solver is written by C++ codes, and Python is used to facilitate pre-processing and post-processing steps and to carry out the stochastic analysis. Here, the polynomial chaos expansions (PCE) method is implemented to quantify the uncertainty of bone mechanics efficiently, and the stiffness matrix is calculated using Hermite polynomials [2]. The new system of equations is numerically solved using Jacobi’s smoother and preconditioned conjugate gradients method. In this study, multigrid and parallelization features are disabled. Young’s modulus and von Mises stress are considered uncertain input and output variables, respectively. Propagated uncertainty through von Mises stress is quantified using statistical analysis of PCE coefficients. The implementation was used to analyse a column with a square cross-section considered a simplified bone to evaluate the correctness of the implementation process. The column meshes using hexahedral elements with 3 elements in width and 10 elements in height. Here, Young’s modulus is homogeneous and normally distributed with mean value and coefficient of variation equal to 1000 MPa and 10%, respectively. Young’s modulus is estimated using a linear combination of two orders of Hermite polynomials. Poisson’s ratio is deterministic and homogeneous equal to 0.3. The bottom face of the column is fixed in the height direction, and the upper face is under displacement-controlled compression (10% of height). Also, rigid body motions are fixed.

Results: Since Young’s modulus is homogeneous and has identical variability through all the elements, it is expected that all stress components will vary identically (i.e. normally distributed with a 10% coefficient of variation). The propagated uncertainty through the predicted von Mises stress is quantified and found to be less than 0.1% different from the analytical result (with respect to mean and variance).

Conclusion: The present evaluation provides an initial demonstration of the correctness of sFE implementation.

Acknowledgements: This study was funded by the EPSRC New Investigator Award (Grant Reference Numbers: EP/V050346/1).

Reference:

[1] La Mattina AA. et al. Ann Biomed Eng 2023; 51:117-124.

[2] Xiu D, Princeton, NJ, USA:Princeton Univ. Press, (2010).

Our research presents a multidisciplinary bio-optical-mechanical approach for describing the three-dimensional (3D) motion of healthy and pathological sperm cells during a free swim in the human female reproductive system. Our novel method is based on a multidisciplinary mechanical-numerical-modeling and experimental-optical approach that combines prior knowledge of the sperm cell's internal structure, and a description of its flagellar beating patterns.

The interaction of human sperm cells with the female reproductive system significantly impacts fertilization. On average, a few dozen to hundreds of sperm cells reach the fertilization site, where only a single sperm cell fertilizes an oocyte. This may be one of the strictest selection processes created by evolution, but its design and role are not fully understood. The prospective correlation between sperm morphology, motility, and uterine interaction should affect the diagnosis and prognosis of male fertility, especially in pathological cases of abnormal sperm morphology and motility. In addition, such analysis may help obtain a more informed sperm selection for in vitro fertilization (IVF), where a male sperm cell fertilizes a female oocyte in a dish. Because the physiological mechanisms of the natural selection of sperm cells in the female body are bypassed in IVF, it is not possible to predict which individual sperm cell would be the one that is most likely to fertilize the oocyte naturally and result in a healthy child, especially in pathological cases. Moreover, the relationship between morphology, motility, cell–surface interaction, and fertilization potential is still not completely clear. Understanding sperm 3D motility while swimming inside the female reproductive tract may lead to improved sperm cell classification and sperm selection for IVF.

We developed a dynamic 3D mechanical finite-element numerical model of sperm cell swimming inside the human female reproductive system and scored the pathological cells according to their chances of reaching the oocyte site. The sperm-cell models, including the full 3D sperm geometry, were preliminary constructed based on experimentally acquired dynamic 3D refractive index profiles of sperm cells swimming freely in vitro as imaged by high-resolution optical diffraction tomography, and then further developed. We numerically stimulated hundreds of cells swimming in a dish and in the female body, with normal and abnormal morphology. We verified that the number of normal sperm cells that succeeded in reaching the fallopian tube sites is greater than the number of abnormal sperm cells. Hence, besides the cell morphology influence, swimming under the female body’s environmental conditions significantly affects the behavior of sperm cells, especially abnormal sperm cells. Moreover, both cell 3D morphology and motility in the complicated geometry of the uterus are significant factors for filtering pathological sperm cells until they reach the oocyte site. Our method may lead to changes in sperm cell classification and evaluation and to new biophysical analysis tools to fill in the gaps from previous studies. Specifically, understanding the migration of normal and pathological human sperm cells inside the female reproductive system can improve the strategies of human sperm cell selection in IVF and fertility evaluation.

Connective tissue is one of the basic tissue types of the body. One of these connective tissues is fascia which is surrounding and interpenetrating skeletal muscle, joints, organs, nerves, and vascular beds. It serves several important functions including transmission of mechanical forces between muscles, tendons, ligaments and bones, protection of organs forming a barrier, regulation of mechanical stress by adsorbing, storing and releasing kinetic energy.

Concerning the mechanical behavior, fascia, is an incompressible, hyperelastic, non-linear and anisotropic material. Anisotropic behavior is given by the spatial orientation of the collagen fibers, which changes along the sheet to ensure a properly response to mechanical demands. As soft tissue, fascia also exhibits viscoelastic properties.

Today, computational simulation is a very powerful tool for the study and analysis of pathologies, treatments and surgeries. To do it correctly, it is necessary an exhaustive characterization and the use of an adequate constitutive model with the ability to predict the behavior of tissues, such as fascia. The present work aims to investigate in depth the mechanical behavior of deep fascia by means of a multidimensional characterization which includes uniaxial (UT), biaxial (BxT) and planar tension (PT) tests. Besides, a constitutive model is proposed to fit tests and obtain material parameters.

Tests were carried out with fascia from sheep. Different shapes were used for testing: dog bone for UT with a central region of interest of 25x5 mm (5:1 aspect ratio), being 25 mm the distance between clamps, rectangular for PT 35x15 mm (width x large), and cruciform shape for BxT with a central region of interest of 15x15 mm. Two directions were defined based on punch orientation with respect to the fibers: longitudinal which has the sample’s axial direction parallel to fibers, and transverse, which has it perpendicular. Testing protocols were defined based on literature and on our own previous tests. For UT and PT three strain levels (2,5%, 5% y 7,5%) are considered performing five cycles in each of them with a strain rate of 10%/min, for BxT case, a 10% strain level is defined undergoing the sample to five cycles and five ratios of one direction to the other, strain rate of 20%/min. Strain measurements and displacements were obtained using GOM Correlate, a digital image correlation software.

To obtain the material parameters by test fitting, a coupled and uncoupled exponential type energy function (SEF) are proposed that consider two perpendicular fiber directions following Stecco 2009.

Results show stress values like to those found in literature that used animal model for testing. For UT, in longitudinal direction, mean stress value was 3,96 MPa, and in transverse 0,6; PT shows for longitudinal and transverse 0,43 and 0,11 MPa respectively; finally, for BxT in case of 1:1 ratio, mean stresses were 3,16 and 1,2 MPa longitudinal to the transverse.

The fitting results show that an uncoupled exponential type function is be able to fit the UT or the equibiaxial experimental data, however fails to predict using the fitting parameters other experiments. On the contrary, the coupled exponential SEF shows good results during the fitting and prediction processes.

Introduction: Microdamage accumulated by cyclic loading or single overloading events contributes to bone fragility through a reduction in stiffness and strength [1]. QCT-based computational modelling cannot incorporate existing in vivo microdamage due to limited resolution. MR imaging on the other hand, is sensitive to pathophysiological changes to adjacent bone marrow that is ‘hidden’ to clinical CT imaging. In the case of repetitive trauma, signal hyperintensity in fluid sensitive sequences is indicative of a stress response where edema, haemorrhage and hyperaemia are present alongside microdamage [2]. Here, we aim to quantify this signal hyperintensity and use it to derive a pre-existing damage variable that represents the underlying tissue damage prior to overloading. We incorporate this variable into a constitutive model to investigate its influence on material and whole bone stiffness and strength.

Materials and Methods: We use the equine athlete as a model for microdamage induced stress fracture where high-speed exercise induces subchondral microdamage. Distal metacarpals (MC3) from n=5 Thoroughbred racehorses were scanned by clinical QCT (0.3 mm voxel size), calibrated to bone mineral density (BMD) and converted to bone volume fraction (BV/TV). MR images (T1w, STIR) were acquired at 3T (0.3 mm voxel size) and registered to the QCT data. Regions of ‘dense’ or ‘sclerotic’ subchondral bone, where microdamage coalesces [3], were segmented from T1w images. Patch-based similarity [4] was used to generate pseudoCT (pCT) images from STIR images. A relative increase in STIR intensity in the dense subchondral bone returned a lower pCT-derived BMD than QCT. This reflects the presence of underlying porosities such as microdamage and increased vasculature [2]. We derived a pre-existing damage variable (D^pex) which affects stiffness and strength and incorporated it into an isotropic, asymmetric BV/TV-dependent elasto-viscoplastic material model (UMAT, Abaqus v6.16) [5,6,7]. Voxel based FEA was used to compress MC3 condyles in silico before (D^pex=0) and after the inclusion of accumulated ‘hidden’ porosity (D^pex>0) to investigate its influence on whole bone mechanical properties.

Results: pCT BMD in dense subchondral bone was lower in all MC3 bones. Incorporating resulted in a median reduction of material stiffness and strength of 20.3% and 20.9% in tension and compression. Inclusion of D^pex reduced whole bone stiffness and strength, and their reduction correlated with D^pex (R²=0.74, R²=0.89). Previous experimental results support our interpretation [8,9].

Discussion: We propose a methodology for incorporating MRI signal hyperintensity into QCT-based FE models to include pre-existing porosities that cannot be detected by clinical CT. Our results illustrate the value of multimodal imaging to potentially capture existing microdamage in vivo. As we use clinical imaging techniques, our results may also aid research in diseases such as osteoarthritis and bone cancer.

References: [1] Seref-Ferlengez, BoneKEy Rep., 4:1, 2015; [2] Taljanovic, Skeletal Radiol., 37:423, 2008; [3] Muir, Bone., 38:342, 2006; [4] Andreasen Med. Phys., 42:1596, 2015; [5] Schwiedrzik & Zysset. BMMB, 12:201, 2013; [6] Schwiedrzik BMMB, 12:1155, 2013; [7] Mirzaali JMBBM, 49:355, 2015. [8] Schwiedrzik Nature Materials 13:740, 2014; [9] Mirzaali Bone 93:196, 2016

Acknowledgements: HBLB VET/CS/027, Siemens Project IPA 42, Leverhulme Trust RPG-2020-215, EPSRC EP/P005756/1.

Neuromusculoskeletal computational models are useful for biomechanical simulations, allowing the analysis of human movement, development of protocols for rehabilitation, physical exercises and training, as well as the development of assistive devices, such as orthoses and prostheses. One of the several types of movements that can be simulated using neuromusculoskeletal models is the sit-to-stand, which is frequently requested during the performance of various activities of daily life of any individual. Despite the importance of simulating such a movement, whether for understanding its execution, analyzing ergonomic aspects, developing assistive devices or identifying possible causes of injury, there is a lack of neuromusculoskeletal computational models dedicated to this type of simulation and that can represent well the biomechanics of a bipedal individual. The OpenSim, an open-source and freely available environment for modeling, simulation and analyzing of the human movement, provides a neuromusculoskeletal model dedicated to squat-to-stand, but such model is not bipedal, which requires the simplification that the individual simulated performs the movement with complete homogeneity in both legs, which does not always happen, especially in cases of hemiparesis, a condition of muscle weakness on only one side of the body due to neuromotor disease (e.g. stroke). Thus, the objective of this work was to develop a bipedal neuromusculoskeletal model capable of representing the human biomechanics of the lower limbs and dedicated to sit-to-stand or squat-to-stand movement simulations. For the development of such a model, we started from the already validated and well accepted model gait10dof18musc freely available by OpenSim. Such model contains 10 degrees of freedom, being able to move in the sagittal plane, and 18 muscles responsible for the execution of the movements. However, such a model is not dedicated to the sit-to-stand movement, being initially designed only for gait. To obtain the model dedicated to the sit-to-stand movement proposed in this work, Hunt-Crossley contact force elements were added between the foot and the ground and kinematic and dynamic constraints were elaborated. In order to test the model obtained for the proposed movement, a predictive simulation was conducted, with the model starting in a sitting position, performing the sit-to-stand and standing stably. Joint movements (hips, knees and ankles) and muscle activations were analyzed. It was verified that the proposed model can perform the sit-to-stand movement well, in a stable way, generating coherent kinematic and dynamic results. The main limitation of the model is related to the foot-to-ground contact, which does not always have enough static for a stable movement, a point that will be improved in later works. In the future, it is intended to improve the model, adding contact force between the hamstrings and a platform capable of representing a seat, improving the foot-to-ground contact ratio and testing it under conditions of hemiparesis.

Breast-cancer and Prostate-cancer are among the most prevalent cancers in women and men, respectively. The World-Health-Organization estimates that about a million deaths occur due to breast and prostate-cancer worldwide due to these cancers each year. Although mostly curable when detected early at the primary site, both of these cancers are incurable when the cancer metastasizes to a distant location in the body, which for these two cancers is eventually bone. There is a scarcity of available human samples and animal models fail due to death preceding bone metastasis; hence a huge unmet need for development of robust in vitro models of bone metastasis. We have developed a novel testbed using a bone mimetic nanoclay scaffold to regenerate human bone followed by sequential seeding prostate and breast-cancer cells obtained from commercial and patient derived cell lines to generate tumors. Extensive analysis of the tumors using gene and protein expression assays and imaging confirm that the testbeds can replicate tumors during mesenchymal-to-epithelial-transition (MET). While many biomarkers exist for evaluation of cancer at the primary site, there are no known bone metastasis markers. Mechanical properties of cells and tissues can capture the complex biological phenomena of adhesion and colonization at the bone site. We also observe via imaging, significant changes to the cytoskeleton quantity and organization within the cells as cancer progresses at bone metastasis. We measured mechanical response of cells using direct-nanoindentation experiments on cancer cells from tumors generated on the testbeds to describe the evolution of cellular properties during cancer progression at the bone metastasis sites. We conducted the nanoindentation experiments under static and dynamic modes to evaluate elastic moduli, hardness, and viscoelastic properties of cancer cells over time. The force-displacement response, elastic moduli, hardness, plastic deformation, viscoelastic properties were captured, with confocal imaging of the cytoskeleton and gene and protein expressions at the same time points. Our results indicate significant reduction in elastic modulus and increased fluid-like behavior of bone metastasized breast-cancer cells (MCF-7) caused by depolymerization and reorganization of F-actin. On the other hand, bone metastasized triple negative cells (MDA-MB-231) showed insignificant changes in elastic modulus and F-actin reorganization over time. We also measured changes to nanomechanical properties of MDA PCa2b(PCa) prostate-cancer cells during the MET and cancer bone metastasis progression over time. The stiffness of PCa cells decreases with metastasis; however, the mechanical plasticity increases during the same time, suggesting that PCa cells become softer on undergoing MET and softer with metastasis progression. In all cases, the imaging and gene, protein expression studies, and computational modeling point towards the depolymerization of actin and reorganization of the cytoskeleton as key factors in the evolution of cell mechanics. In addition, we also subjected the cancer cells to physiologically relevant fluid induced shear stresses that are experimentally enabled using specially designed bioreactors as well as computationally modelled. The role of fluid-derived stresses on migratory characteristics and apoptosis potential of cells is also evaluated. These studies present the use of mechanics-based characterizations as a potential new tool for evaluation of metastasis progression.