Anisotropic inverse bone remodelling algorithm for habitual load case estimation in bone using homogenised finite element models

Associated Publications The theoretical framework for the anisotropic homogenised inverse bone remodelling algorithm is presented in Including Anisotropy in Homogenised Inverse Remodelling Algorithms for Habitual Load Case Estimation in Bone, Gerber and Zysset, Biomech. Model. Mechanobiol., 2026. DOI: 10.1007/s10237-026-02106-1 The theoretical framework for the local anisotropic optimisation scheme is presented in detail in Forward and inverse optimality problems of bone adaptation at the homogenised RVE level, Zysset, Biomech. Model. Mechanobiol., 2026. DOI: 10.1007/s10237-025-02024-8. Project Aim The structure and function of bone are closely related. This project aims to estimate the long-term average loading conditions acting on a bone based on local bone volume fraction and microstructural orientation using time-efficient homogenised finite element (FE) analysis. General Information How to get started The inverse remodelling pipeline requires a running installation of Python and some additional packages. A new conda environment with the required dependencies can be created as follows: conda create --name inv_remodelling python=3.13.2 conda activate inv_remodelling pip install -r requirements.txt Furthermore, the ODB2VTK converter must be installed. Abaqus Standard (Dassault Systèmes) and a Fortran compiler (e.g. Intel oneAPI or GNU Fortran) are required for the FE analysis. The pipeline was tested with Abaqus Standard 2021.HF6. Set path to abaqus command Adjust line 78 in main_inverse_remodelling.py and provide the path to your abaqus command. This is required to run ODB2VTK. How is this repository organised? Code ODB2VTK_Shen2023This folder contains the files for the ODB2VTK. Please visit the official GitHub repository for the newest version. PipelineThis folder contains the main file for the inverse bone remodelling project (main_inverse_remodelling.py) as well as some modules containing custom functions for the individual steps in the remodelling workflow. Further, the UMAT_Maxwell_noPlas.f is used for simulations using a linear elastic material law. Data This folder contains a sample model of a human distal tibia, which serves as a minimal working example for validating the pipeline installation and execution. How can I reproduce the implementation workflow in Gerber et al. 2026? The inverse remodelling pipeline is started by executing main_inverse_remodelling.py. The pipeline will run through the minimal working example saved in the Data folder of this repository. You can rerun the pipeline and verify your installation with the optimisation results provided alongside the sample model. How can I apply the homogenised inverse remodelling scheme to my own projects? You can integrate the optimisation functions from the optimise_stress.py module into your own workflows. The following steps guide you through the inverse bone remodelling process. For detailed information about the individual functions and the input variables, have a look at the function definitions in the corresponding modules. Please note, that the presented code snippets aim to provide you some guidance on a conceptual level and do not provide a working script. For a fully working code sample, have a look at the inverse bone remodelling pipeline for human distal tibia samples main_inverse_remodelling.py. Step 1: Import modules import optimise_stress as opt import post_processing as pp import create_FE_inp as inp_vox Step 2: Prepare input data A dictionary with model information of the following form is required as model input for the optimisation functions. The function inp_vox.parse_input_fractib_MS may provide a helpful starting point on how to import this information from an existing abaqus .inp file. bone = {"Nodes":nodes, # node ID, coord_x, coord_y, coord_z "Elements":elements, # element ID, node IDs for connectivity "bvtv":bvtv, # bone volume fraction "mm":mm, # fabric eigenvalues "orient_line1":orient_line1, # orientation definition from input file (first line) "orient_line2":orient_line2, # orientation definition from input file (second line) "Volumes":vol # element volume } Additionally, a dictionary containing the material constants for the Zysset-Curnier model and the density function f(rho) used in the FE simulations is required. A dictionary containing the material constants used in this project can easily be created by importing the information contained in Data/Sample_Model/1/mat_props.txt using the following method: dict_props = opt.import_input_params(path_to_mat_props.txt) Further, the following variables need to be defined psi_set = 0.037 # value for setpoint complementary free energy density bol_cortex = True # True: cortical elements are included in optimisation scheme; False: cortical elements are excluded from optimisation scheme load_cases = ['fx','fy','fz','mx','my','mz'] # Array listing the load cases for the project bool_lc = [True, True, True, True, True, True] # Boolean array, stating which load cases should be included in the analysis Step 3: Calculate optimal stress state for all elements The optimal stress state for each element based on local bone volume fraction and fabric anisotropy can be calculated as follows: n_elems = bone["Elements"].shape[0] opt_stress = np.zeros([n_elems,6]) for element in range(0,n_elems): # Read out element properties bvtv_element = bone['bvtv'][element,1] mm_element = bone['mm'][element,1:] # Calculate optimal stress opt_stress[element,:] = opt.run_opt_stress(bvtv_element,mm_element,dict_props,psi_set) Step 4: Import simulation results Let's assume you have created six unit load case simulations (fx, fy, fz, mx, my, mz) for your custom model and have run the simulations in Abaqus. As a next step, you need to import the simulation results into Python and store them in a data frame. Ideally, you extract the simulation results from the Abaqus .odb file using ODB2VTK converter. You can then import the results into python and either store them directly in a data frame or save them as parquet files for later use. ## Loop over all load cases for lc in load_cases: # Convert odb to vtu path_odb_res = path_to_odb_file_for_lc pp.odb2vtk_wrapper(path_abaqus,path_odb2vtk,path_odb_res,last_frame) # Determine vtu path path_vtu_file = path_to_vtu_file_created_above # Read stress and strain res_curr = pp.import_vtu(path_vtu_file) cell_data = res_curr.GetCellData() stress_int = cell_data.GetArray('S_Centroid') strain_int = cell_data.GetArray('LE_Centroid') # Store stress and strain stress = [stress_int.GetTuple(el_num) for el_num in range(0,res_curr.GetNumberOfCells())] log_strain = [strain_int.GetTuple(el_num) for el_num in range(0,res_curr.GetNumberOfCells())] # Collect row-wise into flat structure for el_id, (s, e) in enumerate(zip(stress, log_strain)): rows.append({ "sample": sample, "load_case": lc, "element_id": bone_reoriented["Elements"][el_id,0], "volume": volumes[el_id] if el_id < len(volumes) else None, "stress": s, "strain": e }) # Convert to dataframe df_sample = pd.DataFrame(rows) # Save the dataframes df_name = base_dir + folder_FE_models + sample + "/results_" + sample +".parquet" df_sample.to_parquet(df_name,index=False) Step 5: Identify scale factors alpha This step returns the optimised load case as a linear combination of the scaled unit load cases. A, b = opt.OF_derivative_Ab_normalised_lc(bone,df_sample,opt_stress,bol_cortex,bool_lc) alpha = solve(A, b) # solution to A alpha = b Step 6: Calculate objective function values With the function below, you can calculate the objective function values for the optimisation. This allows you the assess optimisation quality. OF_values,OF_values_elements,octant_elements = opt.OF_normalised(bone,df_sample,opt_stress,alpha,bol_cortex) Further considerations The analysis section of main_inverse_remodelling.py and the visualise_results.py module may contain useful functions and methods for the post processing of the optimisation results and for quantitative evaluation of the inverse bone remodelling results for your projects. Detailed Workflow Description The inverse remodelling pipeline is started by executing main_inverse_remodelling.py. At the top of the file, the user can provide project-specific inputs and configure the pipeline. FE Model Generation In this step, the project specific folder structure is created, abaqus input files for 6 canonical load cases are generated, and the FE simulations are run. During the model generation process, the tibia samples are aligned with their anatomical axes in the transversal plane. Import Simulation Results In this step, the abaqus simulation results are imported into Python using the ODB2VTKconverter. Information on sample ID, load case, element ID, element volume, stress and strain are stored in a parquet file for further processing. Load Case Optimisation This step contains the inverse bone remodelling algorithm. First, for each sample the simulation results are imported. In case a St. Venant correction is applied in the optimisation scheme, the elements in the n top- and bottom-most layers are excluded from further analysis. Next, for each element an ideal stress state based on the local bone volume fraction and fabric orientation is determined following the complementary free energy density-based procedure described in Zysset 2026. Then, the scale factors for the canonical load cases are identified by solving a linear equation system on the global scale. A detailed description of the optimisation problem can be found in Gerber and Zysset 2026. The optimisaton results as well as objective function value are stored in a csv output file. Further, objective function values by omitting one load case are calculated.