## 2 Materials and Methods
### Simulation of Bone Formation around the Embedded Interbody Cage  
We assumed that the bone formation around the embedded interbody cage was generated via intramembranous ossification. Intramembranous ossification is an osteogenic pathway and involves osteoclasts, which are cells that break down and resorb bone tissue, and osteoblasts, which are cells involved in bone formation. In many previously published studies in which the bone modelling of the trabecular structure was simulated, bone formation and resorption were calculated in all the elements of the predefined region [15; 16; 20; 22; 23]. In the simultaneous algorithm, bone formation and resorption are calculated simultaneously (Figure 3a). Then, the elements of the formed bone are preserved and the elements of the resorbed bone are removed after the calculation. Thus, an iteration to calculate the density of the elements in the potential region of bone formation and to simulate bone formation and absorption is imposed (Figure 4a). Design space topology optimisation, which describes intermediate structural adaptation progress in the time domain, has been developed and utilised to simulate the complex internal bone architecture of the trabecular bone in the human proximal femur [3; 4; 24]. We adopted the characteristics of design space topology optimisation and developed a sequential bone remodelling algorithm to simulate intramembranous ossification in which layer-by-layer bone formation occurs. In the predefined bone region, the elements adjacent to the superior and inferior vertebrae were imported, and the bone formation and resorption were then calculated as the very first step. After calculation, the elements for the resorbed bone were erased. Then, the elements adjacent to the remaining elements, which were the elements of the formed bone in the previous calculation, were imported, and the bone formation and resorption were calculated again. This process was repeated until new bone formation was no longer generated (Figure 3b). Therefore, in the algorithm, two different iterations were used. One was used to calculate bone formation and absorption, called the internal iteration, and the other was used to import layers adjacent to the vertebrae and generated bone, called the external iteration (Figure 4b). Figure 2: Medical images of the generated bone posterior to the inserted interbody cage after anterior cervical disectomy and fusion (ACDF) surgery. In both the simultaneous and sequential algorithms, the bone formation and resorption were calculated based on the mechanical stimulation of the strain energy density (SED, U) predicted in the FE analysis. The density of the bone, \(\rho\), of the individual elements imported in both the sequential and simultaneous algorithms was calculated using an equation introduced in previously published studies (Equation (2)) [15, 16, 20, 23] as follows:  
\[\frac{d\rho}{dt}=\left\{ \begin{array}{ccc}B(S-(1+s)K)&if&S>(1+s)K\\ 0&if&(1-s)K\leq S\leq(1+s)K\\ B(S-(1-s)K)&if&S<(1-s)K\\ \end{array} \right., \tag{2}\]  
Figure 4: Flow charts of two algorithms. Figure 3: Schematic diagram of bone formation around an embedded interbody cage using two algorithms. where \(dt\), \(B\), \(K\), and \(s\) are the units of time, remodelling rate coefficient, homeostatic equilibrium SED, and threshold range of the lazy zone, respectively. The values for the parameters were set at 30, 1.0 (g/cm\({}^{3}\))\({}^{2}\), 0.0131 J/g, and 0.2443, respectively [25]. The SED per unit mass, \(S\), was calculated as \(S=U/\rho\). The Young's modulus of the individual elements was calculated using Equation (1), and a Young's modulus of the cortical bone of 12,000 MPa was used for the upper boundary. We used 0.1 MPa for the lower boundary, considering the convergence of the FE analysis. Flexion and extension motions were simulated, and the SED, which is the mechanical stimulation used for calculating bone formation and resorption, was predicted using FE analysis.