*2.2. Meshing*

The model's meshing was developed by means of the software ANSYS V.15 (Canonsburg, PA, USA). In summary, the model include 28 cortical bone pieces, 24 trabecular bone pieces, 26 cartilage segments, 4 tendons, 3 ligaments, and the plantar fascia. A trial–error approach was used to optimize the mesh size of each segmen<sup>t</sup> [17]. These authors suggested that the number of inaccurate elements must be less than 5% in all the measured parameters. All simulations and post-processing were developed in Abaqus/CAE 6.14 (Dassault Systèmes, Vélizy-Villacoublay, France) using the available nonlinear geometry solver.

Some of the conditions considered in order to achieve a reasonable mesh size without compromising the calculation time included having a minimum mesh size sufficiently small to fit into the tightest segments, a mesh accuracy of more than 99% of the elements being better than 0.2 mesh quality (Jacobians) and checking that the poor elements were located away from the region of greatest interest (hindfoot bones, metatarsals, PF, and SL) (see Figure 2). The convergence analysis was performed for 265,547 linear tetrahedral elements (C3D4). All parameters exhibited good mesh quality ratios (see Table 1).

**Table 1.** Mesh quality metrics based on Burkhart et al. (2013) recommendations [17].


**Figure 2.** Location of the inaccurate elements, applying the Jacobians as the quality mesh criteria for evaluation.

### *2.3. Tissue Properties*

Two kinds of behavior were considered in this finite element model: linear elastic behavior and hyper-elastic behavior.

Tissues with elastic linear behavior were the cortical bone, trabecular bone, ligaments, and plantar fascia.

Tissues with hyper-elastic behavior were tendons and cartilages.

The numerical values for each of these tissues were as follows:

The material properties (Young's modulus (E) and Poisson's ratio (v)) of the cortical bone, trabecular bone, ligaments, and plantar fascia were assigned in accordance with published data: cortical bone (E = 17,000 MPa, v = 0.3), trabecular bone (E = 700 Mpa, v = 0.3), ligaments (E = 250 Mpa, v = 0.28), and plantar fascia (E = 240 MPa, v = 0.28) [16,20,23].

Tendons and cartilages were modelled as hyper-elastic materials (Ogden model), using the parameters taken from specialized articles [24]. The strain energy density function *U* is:

$$\mathcal{U} = \frac{\mu}{a^2} (\lambda\_1^a + \lambda\_2^a + \lambda\_3^a - 3) + \frac{1}{D} (f - 1)^2,\tag{1}$$

where the initial shear modulus *μ* = 4.4 MPa (cartilage)/33.16 MPa (tendons), the strain hardening exponent *α* = 2 (cartilage)/24.89 (tendons), and the compressibility parameter *D* = 0.45 (cartilage)/0.0001207 (tendons) [16,24,25]. The plantar fascia and spring ligament failures were simulated, applying the isotropic hardening theory that generates a progressive reduction of the tissue's stiffness, resulting in a very flexible material. The initial parameters were a Young's modulus of 240 and a Poisson ratio of 0.3. This strategy allowed us to improve the convergence of the model. Tibialis posterior tendon failure was simulated by removing the traction force of this tendon. Additionally, we considered that this characterization could be more realistic than simply changing the properties of the tissues because it approximates the viscoelastic behavior of these tissues, where stiffness depends on the loading application. The model used for this study maintains the differences of the bone characterization (cortical and trabecular) presented in Cifuentes-De la Portilla, C. et al. [3], where the internal parts of all the foot bones were modeled as trabecular. However, fibula and tibia bones were entirely simulated as cortical because a stress evaluation was not performed on these parts.

### *2.4. Loading and Boundary Conditions*

The FE model was reconstructed from CT images of an unloaded foot. First of all, a standing load position was created (midstance phase) that was used as a reference case to compare against all the pathological cases. In Figure 3, load and boundary conditions are shown. The value of 720 N for the load corresponds to the full weight of an adult of about 70 Kg leaning on one foot. This condition represents a traditional scenario of an AAFD diagnosis assessment. Both loading conditions and boundaries were kept unaltered for all the MCO simulations.

**Figure 3.** Boundary and loading values applied in the model. Loading values correspond to the weight of a person of 70 Kg.

The direction of the load exhibits an inclination of 10 degrees (descending vertical). This load was distributed over the tibia–talus joint (90%) and fibula–talus (10%) [26]. The tendon traction forces were included as reported by Arangio et al. [27]. To simulate the contact with the floor and to avoid the foot structure displacement under loading tests, some nodes located at the lower part of the calcaneus were fixed, while the Z-axis displacement (vertical) of the lower nodes of the first and fifth metatarsals was constrained to 0, using boundary tools available in Abaqus. The nodes remained unaltered for all the simulations performed. To avoid the tendon geometries crossing through the bones, we used the contact surfaces method, using the surfaces of the bones and tendons in contact during simulations.

### *2.5. About the Model Validation*

This study used a model that has been previously validated for other studies related to AAFD [9]. They followed the recommendations of Tao et al. [5], measuring the vertical displacements of some anatomical points: the highest point of the talus (TAL) and of the navicular (NAV), the midpoint of the first cuneiform (CUN), and the highest point of the first metatarsal head (1MT), in two different loading conditions: light loading (minimal contact with the ground) and normal stance loading, using lateral Rx images (sagittal plane) (see Figure 4). The light loading condition was the position before starting simulations.

### *2.6. Model Analysis and Evaluation Criteria*

To quantify the structural deformation of the foot and to evaluate biomechanical stress changes generated by MCO in foot tissues, we simulated the weakness/failure of the plantar fascia (PF), spring ligament (SL) and the tibialis posterior tendon (TPT) in isolation, but also combining these three elements. The weakness/failure was simulated by applying the isotropic hardening theory. This method was implemented using the function "Parameter" in Abaqus, which allows modification of the stiffness of a material, reducing the Young's modulus from its initial value until obtaining a stiffness reduction by about 86%. The stresses in hindfoot bones, forefoot bones, and in all the soft tissues included in the model were calculated.

**Figure 4.** (**Up**) Validation strategy which compares the foot deformation in two different loading values. (**Bottom**) Signs of adult acquired flatfoot deformity achieved with our model.

For measuring stress on tissues, the maximum principal stress (S. Max) was used. This magnitude is closely related to the tensile stress that is generated in foot tissues [28]. Structural deformation was quantified measuring the vertical displacement of the entire structure (in millimeters).
