Finite Element Parametric Design of Hallux Valgus Orthosis Based on Orthogonal Analysis

Objective To design appropriate orthosis for hallux valgus, a difficult foot condition that affects a quarter of the body's bones, we need to clarify the numerical biomechanical features, which have not been established in previous biomechanical studies. Therefore, we constructed a finite element model of the bunion foot to investigate the orthopaedic force compensation mechanism. Methods A patient with moderate hallux valgus was recruited. CT imaging data in DICOM format were extracted for three‐dimensional foot model reconstruction. In conjunction with the need for rapid design of bunion orthosis, a metatarsal force application sizing method based on an orthogonal test design was investigated. The orthogonal test design was used to obtain the hallux valgus angle (HVA) and the inter metatarsal angle (IMA) data for different force combinations. Based on the extreme difference analysis and analysis of variance of the test results, the influence of different force combinations on the bunion angle was quickly determined. Results The results showed that the stress concentration occurred mainly in the first metatarsal bone. The distribution trend was in the medial and lateral middle of the bone and gradually decreased to the dorsal base of the bone body. The greatest stress occurs in the cartilage between the phalanges and metatarsals. In 25 groups of simulation experiments, HVA was reduced from 27.7° to 13°, and IMA was reduced from 12.5° to 7.3°. Conclusion Applying detailed orthopaedic force collocation to the first metatarsal column can effectively restore the mechanics and kinematics of hallux valgus, and provide a reference for the treatment of bunion valgus and the design of orthopaedic devices.


Introduction
A bunion is a foot deformity characterized by a deviation of the bunion and a widening of the first and second intermetatarsal angles. 1 As the medial side of the first metatarsal is often externally compressed by the pressure of footwear, a prominent bump is gradually formed, often accompanied by a cyst, called a bunion, which can cause considerable pain to the patient. 2 Previous biomechanical studies have evaluated metatarsal stress, the valgus angle, and phalangeal displacement after bunion osteotomy, with some studies specifying ligamentous laxity in the foot causing valgus. 3The cause of bunions remains controversial, but overactivity of the metatarsals due to instability of the bunion and bunion muscles near the first metatarsal row is considered the root cause of bunions. 4There are various methods of metatarsal treatment, mainly involving shifting the position of the first metatarsal and changing the length of the first metatarsal to change the pressure distribution in the forefoot and return it to a normal stress state, thus correcting the bunion deformity.
Studies in the year of 2011 have found that the use of bunion orthosis is effective in relieving patients' pain by applying external forces to shift the position of the metatarsal bones to improve the position of the first metatarsal row and prevent pressure. 5Controlling the pronation or eversion by supporting the medial longitudinal arch by orthosis [6][7][8] may slow down the deterioration process and supplement surgical outcomes to reduce recurrence.However, as current orthosis designs rely on subjective perceptions [9][10][11] and do not take into account internal biomechanical influences, they are not as effective as they could be.The reverse-engineered orthosis has a contact surface that fits perfectly to the skin surface, while the internal skeletal anomalies lead to local stress concentrations that can cause redness and pain to the soft tissue surface.Finite elements, because of their ability to model physiological states under different conditions, have been increasingly used in recent years to explore the biomechanical mechanisms of disease and the simulation of treatment modalities.To better match real-life situations, researchers used a finite element foot model reconstructed from normal participants to simulate the potential relationship between altered ligament-to-joint loading on the foot and bunions 12 and constructed a complete foot model for quantitative assessment of the stresses acting on the foot in an upright position. 13The boundary surfaces of the bones and skin of the foot were then analyzed using the finite element analysis software ANSYS.An inverse finite element analysis was performed using an axisymmetric model, and the thickness of the orthopaedic insole was adjusted. 14The potential of finite element simulation for applications in the direction of foot biomechanics is evident.In engineering design, the combination of forces exerted by orthosis is often based on physician experience or subjective patient perception, with large errors detrimental to orthotic outcomes.In the early stages of orthosis design, it is necessary to obtain fast and accurate force application combinations that meet the requirements.Orthogonal tests 15,16 can attain training samples through fewer sample points and reduce the design complexity.This study used finite element simulation to analyze the stress distribution and joint angle variation in the first metatarsal region under the detailed contact surface.It also applied the orthogonal test method to explore the combination and matching of orthotic forces and to investigate the mechanism of orthotic force compensation.
The purpose of this study was to: (i) investigate the biomechanical characteristics of normal and hallux valgus feet under different load conditions; (ii) compare the effects of combination modes of several forces on the load and distribution of the toes and metatarsus; and (iii) provide a reasonable initial value for the optimization iteration of the scheme and a basis for the design of the hallux valgus orthosis by using a simple and low-frequency experimental method.

Finite Element Model Construction
A 22-year-old female patient with a moderate bunion, body mass of 51 kg, height of 160 cm, hallux valgus angle (HVA) of 27.7 , and intermetatarsal angle (IMA) of 12.5 , was selected for CT scanning of the foot.The entire foot below the ankle joint and the initial angle were measured.The CT images were converted into DICOM format and imported into MIMICS V21.0 software (Materialize, Leuven, Belgium), where the toe, tarsal and metatarsal bones, as well as ligaments and muscles, were depicted by the anatomical structure of the human body.The skeletal model was then smoothed and wrapped into a Nubers surface model in Geomagic Wrap 2017 software (3D Systems, Rock Hill, South Carolina, USA), followed by the alignment of the skeletal positions in SOLIDWORKS 2016 software (Dassault Systèmes, Vélizy-Villacoublay, France) and the drawing of the cartilage.1qThe procedure for processing the model is shown in Figure 1.The exact angle value was calculated automatically using the offset angle tool in ANSYS 16.0 software (ANSYS, Canonsburg, Pennsylvania, USA).

Material Properties and Boundary Conditions
The skin material of the foot was simulated with an isotropic and incompressible linear elastic material with a density of 1100 kg/m 3 , 17 an elastic modulus of 0.15 MPa, and a Poisson's ratio of 0.46.The skeletal structures, such as the pelvis, spine, and lower limb bones, were simulated with a rigid material with a density of 1700 kg/m 3 and a Poisson ratio of 0.3.The muscle soft tissue material is very viscoelastic and nonlinear. 18The Mooney-Rivlin model of hyperelasticity in LS_dyna software was used to simulate the strain energy function: where C 10 is the shear modulus, taken as 1.65 KPa; C 01 is the shear modulus, taken as 3.35 KPa; D 1 is the bulk modulus, taken as 3.65 KPa; J el is the elastic volume ratio.
Here, I 1 and I 2 are the first and the second invariants of the Kirschgren equivariant form tensor,  2C.The pressure of the simulated foot model was limited to only moving in the vertical direction to ensure that the angle between the foot and the ground was consistent with the test.Before the simulation began, the foot was placed as close to the bottom plate as possible, without penetrating it, to simulate the slow fall of the foot.Figure 2D shows the simulation results for the human foot sole model.The simulation results of the body pressure distribution cloud for the rigid sole matched well with the test results, with the maximum pressure values located under the heel and forefoot and the body pressure values gradually decreasing outwards along the heel.According to the simulation results, the simulated maximum pressure values were smaller compared to the test values, probably due to the asymmetry of the pressure on the foot when the human body is standing, which leads to an asymmetric distribution of the actual body pressure measured in the test; the presence of the body pressure sensor increases the stiffness of the contact.In summary, the errors in the simulated results for each indicator of the body pressure distribution were small, indicating that the simulated results for the human foot sole model have a high degree of confidence, and the accuracy of the model has been verified.

Finite Element Indicators
To evaluate the effectiveness of the bunion correction, we used two main indicators: the HVA and the IMA.The HVA is the angle between the longitudinal axis of the first metatarsal and the proximal phalanx of the hallux, and the IMA is the angle between the longitudinal axes of the first and second metatarsals.The smaller the HVA and IMA, the better the correction effect.We also calculated the stress distribution and deformation of the foot sole model under different force conditions using ANSYS software.The stress and deformation reflect the biomechanical changes of the foot structure after correction and can provide guidance for optimizing the design of corrective devices.

Orthogonal Test Design and Analysis
Orthogonal design refers to the use of orthogonal tables to rationalize and scientifically analyze multifactor tests, using fewer tests to obtain more accurate results.The orthogonal table uses the number of factors and the level of factors in the test and is arranged according to the principle of orthogonality, which has the characteristics of neat comparability and balanced matching.The orthogonal table can be expressed as L n r m ð Þ, where L represents the orthogonal table code, n is the number of trials, r is the number of levels, and m is the number of factors involved in the trial.Mild bunion foot: 15 < HVA < 20 and/or 9 < IMA < 11 .Moderate bunion foot: 20 < HVA < 40 and/or 11 < IMA < 16 .Severe bunion foot: HVA > 40 , IMA > 20 . 19he currently accepted forces for the correction of bunions are applied from three main positions, as shown in Figure 3: F1 from the distal phalanx of the first bunion to the medial side of the foot, F2 from the area at the metatarsophalangeal joint bunion to the lateral side of the foot, and F3 from the proximal end of the first metatarsal to the medial side of the foot.This paper has a refined contact surface, divided into FS 1 -FS 5 .From the combined simulation analysis of metatarsal force application, the key forces that affect the metatarsal force can be obtained, and the application locations can be subdivided into groups.
The orthotic range of forces under the three-point force orthotic principle is determined by applying moments F1, F2, and F3 to the metatarsophalangeal joint of the foot to provide the initial range of moments for the orthotic optimization of the orthotic forces below.The HVA and IMA were also recorded, and the stress distribution was observed.F1 was loaded from 200 N to 700 N in single increments of 100 N (fix F2 = 300 N, F3 = 400 N); F2 applied pressure values from 100 N to 600 N in single increments of 100 N (fix F1 = 500 N, F3 = 400 N), and F3 applied pressure from 300 N to 800 N in single increments of 100 N (fix F1 = 500 N, F2 = 300 N).
A polar difference analysis is presented in Table 1 to explore the degree of influence of FS 1 -FS 5 on the results from the HVA perspective and the IMA perspective, respectively.By analyzing the magnitude of the extreme differences between the factors, the degree of influence of the different factors on the test results can be obtained.The larger the polar difference, the  greater the influence of the factor on the angle, the most dominant factor.Conversely, it means that the factor has a small effect on the drag coefficient and is a non-dominant factor.

Refinement of Contact Surface Design
Although the foot has a smooth curved surface between the toes, when patients wear commercially available orthosis, the internal bony anomalies, especially at FS 2 , lead to partial pressure concentrations on the skin surface, making them prone to local pressure pain and skin abrasion over prolonged wear.
Based on the conclusions of the metatarsal force application combination simulation analysis, the key forces affecting the metatarsal force are F1 and F2.When F1 force is applied, the division selects the outer end surface of the thumb tip in three parts: S1 is 0-15 mm in length from the tip position, S2 is 15-30 mm from the tip position, and S3 is 30-45 mm from the tip position.
For five different areas, such as FS 1 , FS 2 , FS 3 , FS 4 , and FS 5 , the magnitude of the applied force is used as factors, according to the five levels given in Table 1 using.The experimental design optimizes the extraction conditions, using the proposed level method, to obtain a combined scheme table of 25*5.The table design and the levels of each factor are shown in Table 1, and the maximum stress, maximum displacement, IMA and HVA are calculated at this point.When K i indicates that the level number of a given test factor is Þ , the sum of the results corresponding to k i ¼ K i =5 denotes the arithmetic mean of the resistance coefficient obtained when a factor is taken at level i.
The analysis of the extreme difference results from the HVA perspective, as shown in Table 2. Since the order of polar differences is R FS 3 > R FS 4 > R FS 5 > R FS 1 > R FS 2 , the order of influence for the HVA from primary to secondary is FS 3 , FS 4 , FS 5 , FS 1 , and FS 2 .The two influential factors FS 3 and FS 4 were selected to show the relationship with the HVA, as shown in Figure 4A.The IMA from primary to secondary is FS 1 , FS 2 , FS 4 , FS 5 and FS 3 .The two influential factors FS 1 and FS 2 were selected to show the relationship with the IMA, as shown in Figure 4B.
Combining the relationship between the superior levels of the different factors and the HVA, if there is a positive contribution, e.g., the maximum HVA, the corresponding factor level is selected, that is, the choice of the columns in K i or k i is the corresponding maximum level in each column.However, for bunion metatarsals, the smaller the HVA the better, so picking each K i or k i at the level corresponding to the smallest value is the optimal solution for correcting the HVA.Therefore, FS 1 = 150 N, FS 2 = 100 N, FS 3 = 250 N, FS 4 = 100 N, and FS 5 = 100 N.
As can be seen from Figure 5, the HVA gradually increases for FS 1 , FS 4 , and FS 5 and for FS 2 shows a trend of increasing and then decreasing, but the overall difference is not significant, with FS 3 having the largest fluctuating value and being the optimization focus.Combining the superior level of different factors in relation to the IMA, for bunion orthosis, the smaller the IMA, the better, so picking each K i or K i at the level corresponding to the smallest value is the optimal solution for correcting the IMA.Therefore, FS 1 = 150 N, FS 2 = 100 N, FS 3 = 50 N, FS 4 = 300 N, and FS 5 = 300 N.

Effects of F1, F2, and F3 on Hallux Valgus Angle, Inter Metatarsal Angle, and Stress Distribution
Under the three-point force orthotic principle, the maximum stress occurs in the cartilage in the middle of the proximal phalanx, and its value increases with increasing F1 force, ranging from 83.558 MPa to 292.46 MPa in the range of 200 N-700 N (Figure 6).The pressure distribution also increases with increasing F1 force and tends to increase on both sides.The smaller the HVA, the less desirable the IMA.From 400 N onwards, the HVA stays at around 18 degrees, while the IMA increases from 11 to 15 degrees.We find that the HVA decreases as F1 increases, but continuously increasing the F1 value does not make the HVA continue to decrease; it increases instead.The IMA gradually becomes larger with increasing F1; the recommended range for F1 is 300 N-500 N.
As F2 increased, the stress concentration range of the cartilage in the middle of the proximal phalanx also increased.HVA increased gradually with the increase of pressure in the range of 100 N to 600 N, but HVA exceeded the healthy range when it was above 200 N. IMA gradually decreased with increasing pressure (Figure 7) and reached a healthy range after 400 N. Therefore, the recommended range for F2 is 100 N to 400 N.
The maximum stress was 208.9 MPa in the range of 300 N to 800 N, indicating that F1 was the main factor affecting the maximum stress of the proximal phalangeal bone relative to F3.The pressure distribution of the first metatarsal bone was in the range of 300 N to 800 N, gradually increasing from 8.38 MPa to 12.86 MPa, and progressively increasing with the increase of F3 force.HVA increases gradually with the increase of pressure and changes little after 500 N. Its distribution is concentrated around 22 (Figure 8).IMA decreased gradually with the increase in pressure, but the change was not large.The results showed that compared with F1 and F2, the orthotic effect of F3 was not significant, mainly due to the fixed effect, and excessive force would cause pain in patients.

Orthogonal Test Results
For the optimal solution for the correction of HVA, Group 1: FS1 = 150 N, FS2 = 100 N, FS3 = 250 N, FS4 = 100 N, and FS5 = 100 N; at this time HVA = 12.669 and IMA = 13.689 ; at this time the IMA is large, it is a   , and the stress range was moderate and comfortable (Figure 9).Therefore, in general, the stress is mainly concentrated on the first metatarsal in all four groups.Group 1 has a stress range that is too large, Group 2 is not  suitable for large HVA scenarios, and Group 3 is the average of Group 1 and Group 2, but still not ideal.Group 4 has the best correction effect and the most appropriate stress range.

Summary of Main Findings
In this study, we aimed to find the optimal combination of orthotic forces for correcting mild to moderate bunions using finite element analysis and orthogonal test design.We applied three orthogonal forces (F1, F2, and F3) to the first metatarsophalangeal joint of a bunion foot model under different magnitudes and positions.We evaluated the effects of these forces on HVA, IMA, and stress distribution.We found that F1 and F2 were the most influential factors for HVA and IMA correction, respectively.The optimal combination of forces was

Biomechanical Mechanisms and Factors Involved in Bunion Correction
In the study conducted by Menz et al. 20 on 172 volunteers, it was observed that a higher degree of bunion angle resulted in increased load on the bunion.This phenomenon may be due to the deformation of the bunion, which alters the propulsive function of the first metatarsophalangeal joint, causing an outward deflection of the pressure load on the foot.Our findings are consistent with this, as we have observed higher von Mises stress values in the bunion foot compared to the stress distribution in the first metatarsophalangeal joint region of the normal foot.This explains the inward tilting of the foot and the formation of keratosis in the soft tissue of the medial capsule of the first metatarsophalangeal joint that may occur clinically in bunion patients due to excessive force on the first metatarsal.Reducing the forces on the metatarsals and shifting them towards the bunion is the key to alleviating these problems.
One way to alleviate these problems is to reduce the force on the distal and proximal phalanxes, which can be achieved with certain foot aids.The results indicate that this method can quickly obtain orthotic force combinations that meet requirements, improve the efficiency of analysis, and guide the design of the orthotic shape.Under the principle of three-point force orthosis, the recommended range for F1 is 300-500 N and for F2 is 100-400 N. Therefore, the recommended value for F3 is between 300 and 500 N.It is also known that F1 has a direct, but not linear, effect on the change in the value of the HVA.An appropriate increase in the value of F2 can reduce the IMA while keeping F1 constant.In addition to the correction effect of HVA and IMA, this study also explored the design features and clinical applications of bunion orthosis.We found that bunion orthosis can produce stable corrective force on the bunion through the principle of three-point force.We also found that bunion orthosis can reduce the foot pain and discomfort of patients and improve their quality of life and satisfaction. 5These results are consistent with other studies.Moreover, this study also used a foot stress model in a standing position 13 to analyze the effect of bunion orthosis on foot pressure distribution and found that bunion orthosis can reduce the peak pressure of patients, relieve foot load, and prevent foot complications.

Clinical Implications
Our findings also have important implications for clinical treatment indicators and osteotomy plans for bunions.Rather than solely focusing on reducing HVA, which is often used as a primary outcome measure for bunion correction, indicators should also take into account changes in IMA, which reflect the degree of metatarsal splaying.In the five different directions and positions of the moments considered in this study (FS 1 -FS 5 ), for the HVA, FS 1 , FS 3 , FS 4 , and FS 5 can explain 99.73% of the reasons for the variation and play a dominant influence on the HVA, which can explain why the HVA shows a nonlinear variation when F1 changes, as the counter-directional forces of FS 4 and FS 5 play an important role in limiting lateral flexion.
For the IMA, FS 1 , FS 2 , and FS 4 produced a significant difference in the IMA, with an F-ratio of FS 1 = 13.228 and FS 2 = 15.028,much greater than the F-ratio of FS 4 = 9.536.The force applied to the proximal medial aspect of the phalanx initially reduced the IMA somewhat, but continued increases resulted in an increase in the IMA, increasing the stress on the first metatarsal.These findings provide valuable insights into the underlying mechanisms and factors involved in bunion correction and highlight the need for a comprehensive approach to clinical treatment indicators.These results also have important implications for osteotomy plans, as the aim of metatarsal interception should not only be to facilitate the rotational repositioning of the phalanx but also to achieve external displacement of the first metatarsal to reduce the intermetatarsal angle.Moreover, the study results suggest that orthopaedic aids for bunions should not just push medially on the lateral aspect of the phalanx, as this may not effectively correct bunions and may even increase the intermetatarsal grip and worsen the condition in the long term.

Limitations and Strengths of this Study
This study has some limitations and strengths regarding the use of finite element analysis and orthogonal test design methods for bunion correction: The limitations include simplifying the bone, ligaments, and muscles, unclear experimental conditions, lack of material properties, ligament construction methods, and bunion effects on the first metatarsal and toe joints.The strengths include reducing the number of experiments, refining the contact surface design, performing an extreme difference analysis, validating the simulation results with plantar pressure measurements, and proposing a novel force compensation mechanism.

Conclusion
According to the stress distribution and displacement of the metatarsal bones under various combinations of horizontal forces in this model, the load-bearing capacity of the first metatarsophalangeal joint is reduced due to ligamentous laxity of the joint and high loading pressure.Our findings provide valuable insights into the biomechanical mechanisms and factors involved in bunion correction using orthosis, as well as important implications for clinical treatment indicators and osteotomy plans for bunions.The deformation of the bunion foot can be reduced through the design of biomechanical foot aids to reduce the stress on the metatarsal bones and reduce the patient's pain, providing data to support the design of the orthosis.

FIGURE 1
FIGURE 1 The workflow of a subject-specific musculoskeletal dynamic modeling: (A) Coronal view of CT. (B) Sagittal view of CT. (C) Reconstructed bone tissue.(D) Simulation of the ligament-connecting units.(E) Mesh creation for element analysis.

FIGURE 2
FIGURE 2 Accuracy verification of simulation model: (A) pressure test experiment collection' (B) plantar pressure plate position diagram; (C) test result; and (D) finite element simulation result.

FIGURE 3
FIGURE 3 Diagram of the position of orthotic force application: (A) the method of measuring hallux valgus angle (HVA); (B) direction of orthotic force application; and (C) refining the contact surfaces of the orthotic forces in the F1 and F2 directions.

FIGURE 4 A
FIGURE 4 A plot of HVA angle versus IMA angle.(A) HVA angular relationship diagram; (B) IMA angular relationship diagram.

FIGURE 5 (
FIGURE 5 (A) Horizontal-mean trend graph of FS1 factor for HVA; (B) Horizontal-mean trend graph of FS2 factor for HVA; (C) Horizontal-mean trend graph of FS3 factor for HVA; (D) Horizontal-mean trend graph of FS4 factor for HVA; (E) Horizontal-mean trend graph of FS5 factor for HVA; (F) Horizontal-mean trend graph of FS1 factor for IMA; (G) Horizontal-mean trend graph of FS2 factor for IMA; (H) Horizontal-mean trend graph of FS3 factor for IMA; (I) Horizontal-mean trend graph of FS4 factor for IMA; (J) Horizontal-mean trend graph of FS5 factor for IMA.

FIGURE 6 (
FIGURE 6 (A) HVA and IMA line charts under F1 variation; (B) stress nephogram of the first metatarsal column when F1 is 200N; (C) stress nephogram of the first metatarsal column when F1 is 300N; (D) stress nephogram of the first metatarsal column when F1 is 400N; (E) stress nephogram of the first metatarsal column when F1 is 500N; (F) stress nephogram of the first metatarsal column when F1 is 600N; (G) stress nephogram of the first metatarsal column when F1 is 700N.

FIGURE 7 (
FIGURE 7 (A) HVA and IMA line charts under F2 variation; (B) stress nephogram of the first metatarsal column when F2 is 100N; (C) stress nephogram of the first metatarsal column when F2 is 200N; (D) stress nephogram of the first metatarsal column when F2 is 300N; (E) stress nephogram of the first metatarsal column when F2 is 400N; (F) stress nephogram of the first metatarsal column when F2 is 500N; (G) stress nephogram of the first metatarsal column when F2 is 600N.

FIGURE 8 (
FIGURE 8 (A) HVA and IMA line charts under F3 variation; (B) stress nephogram of the first metatarsal column when F3 is 300N; (C) stress nephogram of the first metatarsal column when F3 is 400N; (D) stress nephogram of the first metatarsal column when F3 is 500N; (E) stress nephogram of the first metatarsal column when F3 is 600N; (F) stress nephogram of the first metatarsal column when F3 is 700N; (G) stress nephogram of the first metatarsal column when F3 is 800N.

FIGURE 9 (
FIGURE 9 (A) The deformation nephogram of foot in Group 1. (B) The strain nephogram of the first metatarsal column in Group 1. (C) The deformation nephogram of foot in Group 2. (D) The strain nephogram of the first metatarsal column in Group 2. (E) The deformation nephogram of foot in Group 3. (F) The strain nephogram of the first metatarsal column in Group 3. (G) The deformation nephogram of foot in Group 4. (H) The strain nephogram of the first metatarsal column in Group 4.
λ 1 is the principal elongation, and γ i is the principal strain.

Table 1
The orthogonal table and the corresponding FE predict the hallux valgus angle (HVA) versus the inter metatarsal angle (IMA) for the first metatarsal column.