Novel Interaction Control in Adolescent Idiopathic Scoliosis Treatment Using a Robotic Brace

Adolescent Idiopathic Scoliosis (AIS) may affect the quality of the patient’s life if it is left untreated. Bracing is prescribed to halt or reduce the curvature progression and avoid surgery. However, the in-brace correction pressure remains unclear, and it is controlled passively by tightening/losing the brace’s strap. Computational modeling has recently attracted researchers’ attention to predict and optimize the AIS bracing treatment. In this paper, a Multi Body-Finite Element (MB-FE) Simscape model and an analytical model of the AIS bracing treatment are created. The MB-FE Simscape model is used to predict the in-brace correction pressure. Furthermore, a Novel Position-based Impedance Control (NPIC) is proposed to control the dynamic interaction between the robotic brace and torso. In this method, the error between the desired and estimated impedance parameters is involved in the controller design to improve the performance of the typical PIC in terms of pose tracking and impedance model tracking. In-vivo data from the literature and numerical simulations are used to validate the MB-FE Simscape model and analytical model. The performance of the proposed controller is verified using numerical simulations in terms of pose tracking and impedance model tracking.


Introduction
Scoliosis is a lateral curvature of the spine, usually combined with a rotation of the vertebrae [1].Adolescent Idiopathic Scoliosis (AIS) is the most common form of scoliosis developed at the ages of 10-18 without any definite causes.AIS occurs in the general population with a wide range of prevalence from 0.93% to 12%; 2 to 3% is the value the most often found in the literature [2].AIS treatment is divided into two main categories: surgical and non-surgical (conservative) treatments.Surgical treatment is typically recommended for patients whose curvatures are greater than 45 while still growing or greater than 50 • when growth has stopped [3].Observation is used for curvatures below 25 • that are still growing or less than 50 • in patients who have completed their growth [3][4][5].Physical therapy [2] and chiropractic medicine [6], classified as alternative treatments, are used for curvatures below 25 • .Bracing, the most common conservative treatment, is used for curvatures between 25 • and 40 [4,5].Several rigid and soft passive braces were developed for scoliosis treatment [7].However, braces have not been integrated with advanced technologies yet, and very few were equipped with smart sensory design and active actuators [8][9][10].The Robotic Spine Exoskeleton (RoSE) with a structure of two Stewart-Gough Platforms (SGPs) was developed by Colombia University [11][12][13].A robotic brace with two SGPs was also developed in France [14].However, they have not been commercialized yet to be used in the clinical treatment.In our previous work, we developed a robotic brace with a structure of three SGPs [7] equipped with 18 linear actuators, 18 position, and 18 force sensors.In addition to the number of SGPs, the significant advantage of our robotic brace is that modeling and control of the physical interaction between our robotic brace and the patient's torso were investigated to improve the safety and compliance of the AIS bracing treatment [7] while the interaction control of the other robotic braces has not been discussed yet [11][12][13][14].
The focus of this study is to propose novel approaches for computational modeling and control of physical interaction in AIS robotic bracing treatment.
First the importance of computational modeling is studied.One of the most critical topics in the AIS bracing treatment is its effectiveness which was widely studied in the literature [15][16][17][18].In [16], 33 AIS patients were followed-up until their skeletal maturity.The results showed that bracing treatment was effective, and 76% of the AIS curvatures were stabilized.In-brace correction pressure directly impacts the AIS curvature correction and the effectiveness of the AIS bracing treatment.Therefore, it is crucial for orthotists to know how much pressure is imposed on the patient's torso by the brace.Some in-vivo studies have been conducted to develop pressure sensory systems to measure the in-brace correction pressure and determine the biomechanical interaction of the typical AIS bracing treatment.A biomechanical evaluation of the AIS treatment using a Milwaukee brace was performed, and the pressure was measured by a commercially available pressure sensor called Dynamic Pressure Monitor (DPM 2000C) [19].The correction pressure exerted on the patient's torso by the Boston brace was measured using the capacitive sensors (PEDAR, Novel, Munich, Germany) [20].Piezoresistive pressure sensors (Tekscan Inc., Boston, MA, USA) were used to measure the correction pressure in the Cheneau brace [21,22].An innovative Smart T-shirt with 100% textile sensors was developed to measure the correction pressure [23].A cost-effective pressure sensing system has recently been built from a piezoresistive polymer placed between two closed-cell foam liners to study the correction pressure exerted by the Cheneau brace on the patient's trunk [24].The pressure reported by the abovementioned in-vivo studies has an extensive range from 1KPa to 124KPa , and the biomechanical interaction between the brace and the patient's trunk is still unclear.Besides, the ethical issue of adjusting the biomechanical performance of the brace using trial and error with human subjects, the long-time needed for adjusting the brace using follow-up and X-ray checks, and concerns of harmful radiation exposure due to repeated radiographic examinations are the main disadvantages of the abovementioned studies.Therefore, computational modeling and numerical simulations have significantly attracted researchers' attention in the last two decades to predict the in-brace correction pressure using virtual models to overcome the limitations of the in-vivo studies.Three main methods have been used to model and simulate the passive bracing treatment of AIS: Multi-Body Modeling (MBM), musculoskeletal modeling, and Finite-Element Modeling (FEM).Although the MB model and musculoskeletal model of the human spine have been widely discussed in the literature [25][26][27], few research works considered the MB model of the AIS bracing treatment.An MB model of the human body was developed to assess efforts along the spine for AIS rehabilitation [28].A musculoskeletal model of the spine and ribcage was developed using AnyBody software to simulate scoliosis [29].A musculoskeletal model of the scoliotic spine was also developed for comparison of ribcage flexibility [30].However, FEM has been widely used to simulate the passive bracing treatment of AIS.FEM analyses were performed to optimize the mechanical structure of the Boston brace [31].The threedimensional corrective force in the AIS bracing treatment was determined using biomechanical FE analysis [32].FEM was also used to predict the biomechanical performance of a textile brace and simulate its effect on a scoliotic spine [33].
The current robotic bracing treatment in [7,[11][12][13][14] are modeled by deriving the analytical model of each SGP independently.The first challenging issue is that the human torso is too complicated to be modeled analytically and numerical methods, e.g., FEM have not been implemented in the robotic bracing treatment of AIS yet.In this paper, a Multi Body-Finite Element (MB-FE) Simscape model and an analytical model of the AIS bracing treatment are created, shown in Fig. 1a.To develop the MB-FE Simscape model, an MB Simscape model of the robotic brace developed in our previous work [7] is first created by importing the 3D SolidWorks model into Simscape.An FE Simscape model of the torso is also created in Simscape using the 'Reduced-Order Flexible Solid' block from the 'Simscape Multibody Library'.The physical interaction between the robotic brace and the human torso is defined using the 'Spatial Contact Force' block available in the 'Simscape Multibody library'.More details of the MB-FE Simscape model are explained in subSection 2.1.
The second issue of current robotic bracing treatment of AIS is that each SGP is modeled and controlled independently [7,[11][12][13][14].To the best of the author's knowledge, the dynamics model of multiple SGPs as a single unit has not been discussed in the literature, while the dynamics model of one SGP has been widely discussed in the literature [34,35].In this paper, the analytical model of the AIS robotic bracing treatment is derived by combining the dynamics formulation of the robotic brace as a single unit and the torso in subSection 2.2.The dynamics formulation of three SGPs as a single unit is derived using the Lagrangian formulation, considering the interaction wrench as an external wrench exerted on the robotic brace.The mass-spring-damper model of the human torso is used to formulate the interaction wrench.The MB-FE Simscape model and the analytical model are validated using the in-vivo stiffness data taken from [13] in subsection 4.1 and Section 4.2, respectively.The desired displacements of the robotic brace required for curvature correction are applied as the primary input to each model and the interaction wrench is obtained from numerical simulations.The displacements of the robotic brace and the interaction wrench are used to estimate the stiffness of each model using the Recursive Least Square (RLS) algorithm.Finally, the estimated stiffness of each model is compared with the in-vivo stiffness data to validate the MB-FE Simscape model and analytical model.Besides, the MB-FE Simscape model is utilized to predict the in-brace correction pressure in sub-Section 4.1 by dividing the interaction wrench by the contact surface area.
In addition to measuring and predicting the in-brace correction pressure, it is also crucial for orthotists to control the in-brace correction pressure/force to increase the efficacy of the bracing treatment.However, the typical braces used in clinics are not equipped with a closed-loop active control system to control the pressure/force and motion exerted by the brace on the patient's torso.The exerted pressure/force is currently controlled passively by regulating the tightness of the brace's strap.Therefore, an active closed-loop control system can improve the efficacy of the current bracing treatment.The motion and force of the robotic braces developed in [11][12][13][14] were controlled using motion control and force control strategies.However, the correction pressure/force and motion of the robotic brace are dynamically dependent.That is why it is vital to control the dynamic relationship between the motion and correction pressure/ force instead of controlling the motion and pressure/force separately.Impedance Control (IC), widely used in robotic applications where the robot has physical interaction with the environment, is a potential solution [36,37].The impedance control is generally divided into two main categories, including Force/torque-based Impedance Control (FIC), also called Impedance Control (IC) in literature, and Position/ velocity-based Impedance Control (PIC) which is also called Admittance Control (AC) in literature.The latter is also classified into the FIC with and without force tracking.The FIC and a model reference adaptive FIC were proposed in our previous work to control the dynamic interaction between the robotic brace and the torso [7].In both FIC and PIC, the dynamic relationship between the motion and the interaction force is modeled as a mechanical impedance, including the desired mass M d , damping C d , and stiffness K d of the inter- action.The FIC regulates the interaction wrench based on the desired impedance parameters according to the motion measurement of the end-effector.A typical force control is then used to provide the regulated interaction wrench.However, the PIC regulates the end-effector motion based on the desired impedance parameters according to the interaction wrench measurement.A typical motion control is then used to provide the regulated motion.The FIC has better performance than PIC when the environment is stiff, while the PIC has better performance than FIC when the environment is soft [38,39].The main limitation of the traditional FIC and PIC is that there are no controller gains to tune the controller performance where the desired impedance parameters M d , C d , and K d are given as constant values.Besides, the actual impedance parameters are not fed back, and the impedance parameters error is not included in the formulation of the control law.In this paper, a Novel Position-based Impedance Control (NPIC), seen in Fig. 1b, is proposed to improve the performance of traditional PIC (AC) in terms of pose tracking and impedance model tracking.The proposed method adds one more feedback loop of the estimated impedance parameters to the typical PIC.The actual impedance parameters are estimated using the RLS estimator and subtracted from the desired impedance parameters to calculate the impedance parameters error.The impedance parameters error multiplied by the proportional gains is used to regulate the interaction wrench.As seen in Fig. 1b, in PIC, the primary input is the measured interaction wrench; however, in NPIC, a correction term obtained from the weighted impedance parameters error is added to the measured interaction wrench to form the primary input of the PIC.In PIC, the virtual pose is computed using the desired impedance model, and a pose control strategy is finally used to achieve the virtual pose.The proportional controller applied to the impedance parameters error gives more freedom to the user to tune the controller performance in terms of pose tracking and impedance model tracking.Section 3 formulates the typical PIC and the proposed controller.In subsection 4.3, numerical simulations are carried out to verify the performance of the proposed control strategy in terms of pose tracking and impedance model tracking using the analytical model derived in subSection 2.2.Although the proposed controller is applied to the robotic bracing treatment of AIS in this paper, it can also be applied to other robotic applications.Table 1 summarizes the analysis of the robotic bracing treatment of AIS in terms of robot design, modeling and control.

Computational Modeling of AIS Bracing Treatment
Computational modeling and numerical simulations can improve and optimize the AIS bracing treatment by modeling and predicting the biomechanical behavior of the AIS bracing treatment.Figure 2 summarizes the main contents of this section.In this section, two models, a Multi Body-Finite Element (MB-FE) Simscape model and an analytical model of the AIS bracing treatment, are developed.Instead of typical braces being used in the AIS treatment, the robotic brace developed in our previous work [7] is modeled to mimic the biomechanical behavior of typical braces.In the MB-FE Simscape model, an MB Simscape model of the robotic brace is created by importing the CAD model into Simscape.
An FE Simscape model of the human torso is also created using the 'Reduced-order Flexible Solid' block from the 'Simscape Multibody Library'.Finally, the 'Spatial Contact Force' block from the 'Simscape Multibody library' is used to model the physical interaction between the robotic brace and the torso.In addition, the analytical model is created by deriving the dynamics formulation of the robotic brace with three SGPs as a single unit using the Lagrangian formulation, assuming that the interaction wrench obtained from the mass-spring-damper model of the torso is added as an external wrench to the dynamics model of the robotic brace.The MB-FE Simscape model and analytical model are validated in subSections 4.1 and 4.2, respectively.The MB-FE Simscape model is used to predict the in-brace correction pressure in subSection 4.1.In subsection 4.3, the analytical model is used to simulate the interaction control of the AIS bracing treatment and verify the proposed interaction control strategy which will be formulated in subSection 3.2.

MB-FE Simscape Model
Figure 3 shows the schematic diagram and the visualized form of the MB-FE Simscape model.As seen in Fig. 3a, three main steps are followed to create the MB-FE Simscape model: robotic brace modeling, human torso modeling, and physical interaction definition.The first step is to create an MB Simscape model of the robotic brace using its SolidWorks model.The robotic brace was described in details in [7].create the FE Simscape model of the torso.The 3D scan of a scoliosis patient used in our previous work [7] is saved as a '.stl ' file, and the '.stl' file is imported into Simscape using this block.FE analysis in MATLAB is used to derive the torso's reduced-order model dataset, including the torso's stiffness, damping, and mass matrices.The following procedure is followed to perform the FE analysis and generate the reduced-order model dataset: 1) the 3D scan of the patient's torso is used in this FE analysis.2) the material properties adopted from the literature are Young's modulus E = 10000Pa [40], Poisson's ratio = 0.394 , and the mass density = 1250kg∕m 3 [41]. 3)the locations of the interface frames are defined as [0 0 0; -7.41 266.91 -27.33; -8.63 507.4 -50.17] mm.The interface frames determine the connection points used to define the contact surfaces of the torso model.4) the finite element mesh is generated using MATLAB's command 'generateMesh'.5) the multipoint constraints for the interface frames are defined using MAT-LAB's command 'structuralBC' to preserve the six degrees of freedom at each frame.6) the reduced-order dataset is obtained using MATLAB's command 'reduce'.The visualized form of the FE Simscape model of the torso is seen in the middle box of Fig. 3b.
As seen in Fig. 3a, the third step is to define the physical interaction between the MB Simscape model of the robotic brace and the FE Simscape model of the torso.The 'Spatial Contact Force' block is used to define the physical interaction.This block applies normal and frictional contact force between the torso and the moving rings.However, it is impossible to directly connect the 'File Solid' blocks used for modeling the torso and the moving rings to the 'Spatial Contact Force' block.Therefore, contact proxies are used to simulate the contact surfaces.18 'Solid Bricks' are connected as contact proxies to the reduced-order flexible model of the torso using 18 'Rigid Transform' blocks to represent the torso's contact surfaces.18 'Solid Brick' blocks are also rigidly connected to the moving rings in the MB Simscape model to represent the contact surfaces of the moving rings.The proxies are directly connected to the 'Spatial Contact Force' blocks.Note that the physical interaction is defined between the entire geometry of the proxies by ticking the 'Entire Geometry' checkbox in the dialog box of each 'Solid Brick' block of the proxies.Some of the torso's proxies are seen in the middle box of Fig. 3b.The right-hand side box in Fig. 3b also represents the visualized form of the MB-FE Simscape model.

Analytical Model
In this subsection, the analytical model of the AIS bracing treatment is derived by formulating the dynamics model of three SGPs as a single unit, while the interaction wrench F e is defined as an external wrench exerted on the moving rings.The schematic diagram of the robotic brace is shown in Fig. 4a in which the geometrical parameters of the first SGP are seen.The architecture of the robotic brace is described first.It consists of three SGPs (four rings) connected in series.The numbering of the rings, SGPs, and actuators are the same as what is seen in the right-hand side of Fig. 3b.It has 18 active DOFs controlled by 18 linear actuators.The architecture of each SGP follows the kinematic structure of a 6-6 SGP, where all the limbs have identical kinematic chains of SPS.The global frame {A} is attached to the fixed ring at the origin O A and frames {B}, {C} , and {D} are attached to the first, second, and third moving rings at the origins O B , O C , and O D , respectively.The attachment points of the first, second, and third SGPs are described by A i1 , B i1 , A i2 , B i2 , and A i3 , B i3 , respec- tively.The position of the attachment points of the first, second, and third SGPs with respect to (w.r.t.) their local frames are denoted by respectively.The length of the i th limb of the first, second, and third SGPs are represented by l i1 , l i2 , and l i3 , and the unit vectors ŝi1 ,ŝ i2 , and ŝi3 are used to describe the direction of the corresponding limbs.The position of the origins O B , O C , and O D w.r.t. the fixed frame {A} are described by the position vectors A P 10 = [p 1x p 1y p 1z ] T , A P 20 = [p 2x p 2y p 2z ] T , and A P 30 = [p 3x p 3y p 3z ] T .The orientation of the moving rings w.r.t. the frame {A} is also described using screw coordi- nates representation as T represents the velocity of the first, second, and third moving rings. 1 , 2 , and 3 denote the angular velocity of the first, second, and third moving rings, respectively.The linear velocity of the first, second, and third moving rings are defined as A v 10 , A v 20 , A v 30 , respectively.
Figure 4c indicates the main steps followed to derive the analytical model.The method proposed in our previous work [7] to formulate the dynamics model of one SGP is extended to derive the dynamics formulation of three SGPs as a single unit in the presence of physical interaction.First, the dynamics formulation of the limbs is derived using the Lagrangian formulation.Second, the dynamics formulation of the moving rings is derived using the Lagrangian formulation, considering that the interaction wrench F e is assumed as an external wrench exerted on the moving rings.Finally, the dynamics formulation of the moving rings and limbs are added together to form the analytical model of the AIS robotic bracing treatment as follows: where M( ) , C(, χ ) , and G( ) denote the 18 * 18 mass matrix, the 18 * 18 Coriolis and centrifugal matrix, and the 18 * 1 gravity vector of the robotic brace.J and denote the 18 * 18 Jacobian matrix and 18 * 1 actuator forces vector.F d and F are 18 * 1 vectors denoting the disturbance wrench and the projection of the actuator forces.
The dynamics matrices of the robotic brace is obtained as follows in which M p , C p , and G p denote the 18 * 18 mass matrix, the 18 * 18 Coriolis and centrifugal matrix, and the 18 * 1 gravity vector of the moving rings.M li , C li , and G li represent the 18 * 18 mass matrix, the 18 * 18 Coriolis and centrifugal matrix, and the 18 * 1 gravity vector of the i th limb of the first, second, and third SGPs for i = 1, 2, … , 6 .The dynam- ics matrices of the limbs are as follows in which K li and P li denote the kinetic and potential energy of the limbs.M ′ is obtained from kinetic energy as follows (5) where m pi and M pi denote the mass and mass matrix of the i th moving ring.K pi and K p represent the kinetic energy of the i th moving ring and robotic brace.P pi and P p represent the potential energy of the i th moving rings and robotic brace.I pi indicates the inertia of the i th moving ring.The dynamics modeling of the human torso is too complicated.Since non-linear control is not our focus in this paper, we made an assumption in which a linear massspring-damper model of the torso is used to compute the interaction wrench.18 mass-spring-damper systems are used to model the thoracic segment of the torso at the vertebrae T11, T7, and T4 levels, seen in Fig. 4b.Therefore, the interaction wrench is formulated using the mass-springdamper model of the torso as follows in which M T11 , C T11 , and K T11 are 6 * 6 diagonal matrices that denote the torso's mass, damping, and stiffness matrices at the vertebrae T11 level.The 6 * 6 diagonal matrices M T7 , C T7 , and K T7 represent the mass, damping, and stiff- ness matrices of torso at the vertebrae T7 level.The mass, damping, and stiffness matrices of the torso at the vertebrae T4 level are represented by 6 * 6 diagonal matrices M T4 , C T4 , and K T4 .e1 also indicates the initial position of ( 21) the interaction points.Besides, the interaction wrench is T F e,T7 T F e,T4 T ] T in which F e,T11 , F e,T7 , and F e,T4 denote the interaction wrench vector at T11, T7, and T4 levels, respectively.F e,T11 is computed as the resultant wrench of six individual interaction wrenches of each mass spring damper model ( F e1 , …, F e6 ), seen in Fig. 4b.F e,T7 , and F e,T4 are computed in the same way as F e,T11 .

Interaction Control of AIS Bracing Treatment
As mentioned in Section 1, an active control strategy is needed to provide the required correction pressure/force for AIS curvature correction.The correction pressure/ force and motion of the torso are dynamically dependent.Therefore, the Force-based Impedance Control (FIC) and a model reference adaptive FIC were proposed in our previous work to control the dynamic relationship between the correction pressure/force and motion instead of controlling the correction pressure/force and motion separately [7].In this section, a Novel Position-based Impedance Control (NPIC) is proposed to improve the performance of the typical PIC in terms of pose tracking and impedance model tracking.The typical PIC is introduced in subSection 3.1, and the NPIC is formulated in subSection 3.2.

Position-Based Impedance Control (PIC)
The schematic diagram of the PIC without force tracking is shown in Fig. 5a.The mass-spring-damper system is used to model the mechanical impedance and represent the dynamic relationship between the robotic brace and the environment.Therefore, the desired impedance model is defined as follows in which F d e and e = d − represent the desired interac- tion wrench and pose error, respectively.d also denotes the desired pose.In PIC, first, the virtual pose of the robot's end-effector v required for providing the desired impedance model in Eq. ( 27) is computed by rewriting the desired impedance model in Eq. ( 27) as follows, while the desired interaction wrench F d e in Eq. ( 27) is replaced with the measured interaction wrench F m e and the desired pose d is given.
Second, typical pose control, e.g., Inverse Dynamics Control (IDC), is used to guarantee that the virtual pose of the robot's end-effector v is achieved as follows ( 27) Finally, F τ is multiplied by the inverse transpose of the Jaco- bian matrix of the robotic brace J −T to obtain the required actuator forces .To add the force tracking ability to the typical PIC, the desired interaction wrench F d e in Eq. ( 27) is replaced with F Ref e − F m e in which F Ref e denotes the reference interaction wrench.Therefore, the virtual pose v is computed by rewriting Eq. ( 27) as follows

Novel Position-Based Impedance Control (NPIC)
The main challenge of the PIC is that there are no control gains to tune the PIC performance in terms of pose tracking and impedance model tracking, while the desired impedance parameters M d , C d , and K d , are defined as constant values.Note that K p_IDC and K d_IDC in the inner loop of the PIC are only used to achieve the virtual pose.The question is how to tune the virtual pose to achieve the desired performance in terms of desired pose tracking and impedance model tracking.The solution to this problem is coming up by considering the basic principles of the typical position control and force control strategies.In typical position control and force control, the position error and force error signals are used to formulate the control law, respectively.However, in traditional impedance control, the impedance parameters error is not involved in formulating the control law.In other words, there is no feedback loop of the estimated impedance parameters.The main concept of the proposed controller is to add a new feedback loop of the estimated impedance parameters to the typical PIC and involve the impedance parameters error in the control law formulation.
The schematic diagram of the NPIC is shown in Fig. 5b.First, the actual impedance model is formulated as follows in which Md , Ĉd , and Kd represent the actual mass, damping, and stiffness of the interaction.The Recursive Least Square (RLS) estimator is used to estimate the actual impedance parameters while the measured interaction wrench and pose error signal are applied as inputs to the estimator.The actual impedance model is rewritten as follows in which  T (t) = ë(t) ̇e(t) e(t) and ̂ im = Md Ĉd Kd T denote the regressor and the estimated impedance (30) parameters vector.The detailed formulation of the RLS method is as follows in which L k , P , and denote the correction weight, covari- ance, and forgetting factor of the estimator.Second, proportional controllers with the gains K 1 , K 2 , and K 3 are applied to the impedance parameters error to obtain the weighted impedance parameters error.K 1 , K 2 , and K 3 are 18*18 diago- nal matrices as follows in which K 1_T11 , K 1_T7 , and K 1_T4 represent the proportional gains for the first, second, and third moving rings multiplied by the mass error M d − Md .K 2_T11 , K 2_T7 , and K 2_T4 also denote the proportional gains for the first, second, and third moving rings multiplied by the damping error C d − Ĉd .The proportional gains for the first, second, and third moving rings multiplied by the stiffness error K d − Kd are repre- sented by K 3_T11 , K 3_T7 , and K 3_T4 .Third, the impedance model error is computed by multiplying the weighted impedance parameters error by the pose error e and its derivatives.The computed signal is added to the measured interaction wrench F m e to form the desired interaction wrench F d e as follows Finally, the desired interaction wrench F d e obtained in Eq. ( 39) is applied as the primary input to the typical PIC, seen in Fig. 5b.Note that the measured interaction wrench F m e is the primary input to the typical PIC, seen in Fig. 5a.Therefore, the virtual desired pose x v in NPIC is computed using Eq. ( 28) where F m e in Eq. ( 28) is replaced with F d e in Eq. (39).The IDC's control law in Eq. ( 29) is then used to control the pose of the robotic brace such that the virtual desired pose x v is achieved.Although the NPIC is formu- lated in 18 DOFs for the AIS bracing treatment, it can also be applied to other robotic applications.
In comparison to FIC and PIC, the NPIC regulates both interaction wrench and the virtual pose of the end-effector according to the impedance parameters error and measured interaction wrench, respectively, while the FIC only regulates the interaction wrench according to the measured endeffector pose and the PIC only regulates the end-effector pose according to the measured interaction wrench.Besides, the NPIC gives more freedom to the user to improve the performance of the typical PIC in terms of pose tracking and impedance model tracking by tuning the proportional gains K 1 , K 2 , and K 3 applied to the impedance model error.

Numerical Simulations
In this section, the MB-FE Simscape model and analytical model of the AIS bracing treatment are validated using numerical simulations and the in-vivo data taken from [13].The in-brace correction pressure is also predicted using the MB-FE Simscape model.Besides, the proposed controller is verified using numerical simulations by applying the NPIC to the analytical model.Figure 6 summarizes the main contents of this section.First, the MB-FE Simscape model of the AIS bracing treatment is validated by comparing the estimated stiffness coefficients obtained from numerical simulations in Simscape with the in-vivo stiffness data from [13].The inverse kinematics problem is solved to compute the actuator displacements in the joint space based on the given desired displacements of the moving rings in the task space because the spatial motion of the moving rings is created by the active prismatic joints in the joint space.Then, the actuator displacements are applied as the inputs to the MB-FE Simscape model, and the output interaction wrench is obtained.The interaction wrench and the desired displacements of the moving rings are imported into MATLAB's Curve fitting Toolbox (CFT), and the slope of the linear regressor fitted to the measurements using the Recursive Least Square (RLS) algorithm is considered as the estimated stiffness of the thoracic region at T11, T7, and T4 levels.The in-brace correction pressure is also predicted by dividing the interaction wrench by the contact surface area.Second, the analytical model is validated by comparing the estimated stiffness of the analytical model obtained from MATLAB simulations with the in-vivo stiffness data from [13].The desired displacements of the moving rings and the actuator forces obtained from the previous Simscape simulation are applied as the inputs to the analytical model.The same procedure as the previous Simscape simulation is followed to estimate the stiffness of the analytical model.Third, the performance of the NPIC is compared with PIC in terms of position tracking and impedance model tracking to verify the NPIC using the numerical simulations in MATLAB.The geometric and inertial data of the first, second, and third SGPs used in this section are described in Table 2.

Validate the MB-FE Simscape Model and in-Brace Correction Pressure Prediction
Numerical simulations in Simscape are carried out to validate the MB-FE Simscape model by comparing the simulation results with the in-vivo data.Then the in-brace correction pressure is predicted using the MB-FE Simscape model, assuming that the desired displacements of the torso required for moving the spine to the desired posture are given.The in-vivo stiffness data from [13] is used to validate the MB-FE Simscape model.In [13], eight healthy males and two individuals with spine deformities participated in the study to take a series of force-displacement measurements.The Robotic Spine Exoskeleton (RoSE) was used to apply six unidirectional displacements in each DOF on the thoracic and lumbar segments of the participants, including −15, −10, −5, 5, 10, and 15mm for translation and  The position/orientation and force/moment were measured, and the stiffness matrices were derived.In our simulation, −12mm displacement along the x direction in the coronal plane is considered as the desired displacement of each moving ring Δx d .Note that the desired trajectory is assumed to be a cubic polynomial [35].Since the spatial motion of the moving rings is created using the linear actuators modeled by prismatic joints in the MB-FE Simscape model, the inverse kinematics problem of the robotic brace is first solved to compute the desired displacements of the prismatic joints.The desired displacements of the prismatic joints are applied as the inputs to the MB-FE Simscape model, and the interaction force is obtained.The stiffness coefficients of the MB-FE Simscape model in the x direction are estimated using the RLS estimator in MATLAB's CFT as follows: The interaction force and torso displacement measurements in the x direction are imported into MATLAB's CFT, and a linear model is fitted to the data.The linear model's slope is considered the estimated stiffness of the MB-FE Simscape The estimated stiffness of the MB-FE Simscape model are finally compared with the in-vivo stiffness data from [13] to validate the MB-FE Simscape model.Note that the human subject in [13] is in the sitting posture, therefore, the hip and lumbar segments of the torso model is fixated on the fixed ring of the robotic brace to simulate the sitting posture of the torso model.
The comparison between the stiffness of the MB-FE Simscape model and the in-vivo data is shown in Table 3.The stiffness coefficients are estimated at three levels of the thoracic region, including T4, T7, and T11 vertebrae.In [13], the stiffness of the thoracic region in the x direction is 2758.34± 361.22N∕m .Note that the stiffness of the tho- racic region obtained in [13] has a large standard deviation of 361.22N∕m .One reason is that the stiffness matrices in [13] were computed by taking the average of the stiffness matrices of eight individuals in which the stiffness of the torso for each person is entirely different from the other participants.According to Table 3, the stiffness coefficients of the MB-FE Simscape model in the x direction are 2933 , 2714 , and 2730N∕m for T11, T7, and T4 levels, respectively.Besides, the errors between the stiffness of the MB-FE Simscape model and in-vivo data are 6.33% , 1.61% , and 1.02% for T11, T7, and T4 levels, respectively.Therefore, it is concluded that the stiffness properties of the MB-FE Simscape model are in the same range as the in-vivo data from [13].Note that since the main focus of this paper is scoliosis treatment, and the spine curvature of scoliosis occurs mainly in the coronal plane; the stiffness coefficients in the coronal plane ( x direction) are studied in these simulations.
The MB-FE Simscape model is also used to predict the in-brace correction pressure required for AIS curvature correction.The mean interaction force obtained from the numerical simulations in Simscape is divided by the rectangular contact surface area of the MB-FE Simscape model ( 6cm 2 ) to compute the predicted in-brace correction pres- sure.The correction pressure of the thoracic region at the T11, T7, and T4 levels are calculated as 31.6,91.6 , and 62.5Kpa , respectively.Therefore, the mean value for the pre- dicted pressure exerted on the thoracic region is 63.86Kpa .The reported pressure range in [19][20][21][22][23] is from 1Kpa to 112Kpa .The pressure data reported in the most recently published paper [24] ranges from 14Kpa to 112Kpa .It is concluded that the predicted pressure obtained from the MB-FE Simscape model is in the same range as the pressure data reported in the literature.Note that the reasons behind the extensive range of the pressure reported in the literature are that different braces, including Milwaukee, Boston, and Cheneau, different patients, and different sensing systems are used in each research work for measuring the pressure.Besides, the type of activities, including standing, walking, supine, and laying, also affects the pressure ranges.

Validate the Analytical Model
The analytical model presented in Eq. ( 1) is validated by comparing the estimated stiffness of the analytical model obtained from MATLAB simulations and the in-vivo stiffness data from [13].It means that the dynamics matrices of the robotic brace M(), C(, χ ), and G( ) derived using the Lagrangian formulation and the assumption of adding the interaction wrench as an external wrench to the dynamics formulation of the robotic brace to form the analytical model are validated.Therefore, the inputs and outputs for the validation process must be chosen from three main quantities existed in Eq. ( 1), including the motion variables ( , χ , and χ ), the projection of the actuator forces F = J T , and the interaction wrench F e .In this simulation, the motion variables and the projection of the actuator forces are considered as the inputs and the interaction wrench is defined as the output of the analytical model.The desired motion of the analytical model is considered a cubic trajectory with a −12mm displacement along the x direction in the coronal plane for each moving ring.The actuator forces required for creating such a motion are taken from the numerical simulations in subSection 4.1.Therefore, the interaction wrench is computed as follows Now the Curve Fitting Toolbox in MATLAB is used to fit a linear regression model on the interaction wrench F e and displacement measurements Δx d using the RLS algo- rithm.The slope of the linear regression model is considered the estimated stiffness of the analytical model, and it is compared with the in-vivo stiffness data from [13].The comparison between the estimated stiffness of the analytical model and the in-vivo stiffness data of the thoracic region is shown in Table 4.The estimated stiffness of the analytical model at the T11, T7, and T4 levels in the x direction are 2829N∕m, 2860N∕m, and 2645N∕m while the stiffness of the thoracic region in the x direction reported in [13] is 2758.34+ 361.22N∕m .The errors between the analytical model's stiffness at the T11, T7, and T4 levels and the in-vivo stiffness data are 2.56%, 3.69%, and 4.11% .It is concluded that the stiffness coefficients of the analytical model have an acceptable range with respect to the in-vivo stiffness data.
One limitation of the MB-FE Simscape model is that the stiffness and damping properties of the 'Spatial Contact Force' block affect the interaction force in the MB-FE Simscape model.This also affects the interaction force obtained from the analytical model because the actuator forces for the analytical model are taken from the MB-FE Simscape model.Although some differences between the interaction force obtained from both models (MB-FE Simscape model and analytical model) and the interaction force reported in ( 40) [13] may exist, the predicted interaction pressure is still in the same range as the other literature ( 1Kp to 112Kpa [19][20][21][22][23][24]). Note also that one reason behind this difference is that only one moving ring is used for manipulating the thoracic segment in [13], while we used three moving rings for the thoracic part and all the rings are moving at the same time in the simulations.Besides, the interaction force obtained from models can be reduced (increased) by reducing (increasing) the stiffness and damping coefficients of the 'Spatial Contact Force' block to be in the same range as [13].
It is suggested to use Digital Twin (DT) Technology and Machine Learning (ML)-based parameter identification methods to create a Digital Twin of the AIS bracing treatment and optimize the stiffness and damping parameters of the 'Spatial Contact Force' block while the real-time force-displacement measurements are collected using realtime experiments for the parameter identification.Therefore, the stiffness and damping parameters of the 'Spatial Contact Force' block are tuned such that the error between the interaction force obtained from the DT and the real-time experiments is zero.Now we are also working on creating a DT of the AIS bracing treatment using ML-based parameter identification methods and real-time experiments to update the stiffness and damping parameters of the 'Spatial Contact Force' block.

Verify the proposed controller
Numerical simulations in MATLAB are carried out to verify the performance of the NPIC in terms of position tracking and impedance model tracking.To this goal, the NPIC and typical PIC are applied to the analytical model, and the performance of the NPIC is compared with the PIC in terms of position tracking and impedance model tracking.It is assumed that the physical interaction occurs only in the x direction, and the interaction wrench in the other DOFs is considered zero.Therefore, the NPIC and PIC are applied to each moving ring only in the x direction, and typical Inverse Dynamics Control (IDC) is used to control the pose of each moving ring in the other DOFs.The desired displacements of the moving rings in the global x direction are assumed −12mm and for the other DOFs 0mm .The estimated impedance parameters obtained using the RLS estimator from the NPIC and PIC are compared with the desired impedance parameters to investigate if the desired impedance model is reached.The simulated interaction force in the x direction is created using Eq.(26).The stiffness coefficients of the thoracic segment of the torso model in the x direction at three contact surfaces, includ- ing the vertebrae T4, T7, and T11 levels, are considered K T11_x = K T7_x_ = K T4_x = 2758.34N∕mtaken from [13].The damping coefficients in the x direction are also assumed The desired impedance parameters of the interaction in the x direction for each mov- ing ring are defined as The proportional and derivative gains of the IDC ( K p_IDC and K d_IDC ) have the diagonal elements of 1000 and 10 for the translational DOFs, respectively.The proportional and derivative gains of the IDC K p_IDC and K d_IDC ) for the rotational DOFs are defined as 10 and 1 , respectively.
The design criteria is to reduce the mean percentage error of position tracking and impedance parameters tracking to be smaller than 4% ( 0.48mm ) and 1% , respectively.To achieve this goal, two different sets of gains for the NPIC are used.First, it is assumed that the proportional gains of the NPIC for three moving rings are 7 represents the desired position and the position of three moving rings in the x direction obtained from NPIC and PIC.In NPIC, the steady state error of the position tracking is too close to zero, while the PIC has a large position tracking error.In contrast to typical PIC, the proposed controller gives the robotic brace the capability of tracking the desired trajectory in the presence of the physical interaction between the robotic brace and torso.Besides, the estimated impedance parameters obtained from NPIC and PIC for three moving rings along with the desired impedance parameters are shown in Figs. 8,  9, and 10.It is seen that the desired stiffness and damping behavior are reached, although the estimated mass for both NPIC and PIC is not close to the desired mass, and the Fig. 7 Position tracking of the moving rings in the x direction for the first set of gains (Desired-NPIC-PIC): the position error of NPIC is much smaller than that of PIC Fig. 8 Comparing the estimated impedance parameters obtained from NPIC and PIC with the desired values for the first set of gains (first SGP): the tracking error of stiffness and damping in NPIC is less than those of PIC, while mass tracking error in PIC is smaller than that of NPIC estimated mass obtained from the NPIC has a bigger estimation error than the estimated mass obtained from PIC.However, the proportional gains of the proposed controller K 1 , K 2 , and K 3 can be tuned such that the mass estimation error goes to zero and the desired impedance model is achieved.The second set of the proportional gains of the NPIC is assumed as  5 indicate the mean percentage error of the position in the x direction and impedance parameters for two set of gains, while the rows represent the mean percentage Fig. 9 Comparing the estimated impedance parameters obtained from NPIC and PIC with the desired values for the first set of gains (second SGP): NPIC provides better performance than PIC in terms of desired stiffness and damping tracking, while PIC has better performance than NPIC in terms of desired mass tracking Fig. 10 Comparing the estimated impedance parameters obtained from NPIC and PIC with the desired values for the first set of gains (third SGP): in contrast to PIC, the desired stiffness and damping is followed well using NPIC, but the NPIC performance needs to be tuned in terms of mass tracking error of the first, second, and third SGPs for NPIC and PIC.The first set of gains is discussed first in terms of position tracking and impedance model tracking.The main goal of control strategies proposed in Section 3 (PIC and NPIC) is to reduce the position tracking error and provide the desired mechanical impedance between the robotic brace and torso (follow the desired impedance model).As seen in the first column of Table 5, the mean percentage error of position for the first, second, and third SGPs obtained from NPIC (3.2797%, 1.2950%, 0.8755%) is less than those of PIC (16.9347%, 11.4851%, 10.5578%).The mean percentage error of the damping and stiffness parameters for the first, second, and third SGPs obtained from NPIC is also less than those of PIC, seen in the fifth and seventh column of Table 5.But the mean percentage error of the mass for the first, second, and third SGPs obtained from NPIC is bigger than those of PIC, seen in the third column.In brief, the proposed controller NPIC improves the performance of PIC in terms of position tracking and desired stiffness and damping tracking, while the NPIC cannot follow the desired mass well.
The second set of gains is then discussed in terms of position tracking and impedance model tracking.The mean percentage error of the mass for the first, second, and third SGPs obtained from NPIC is much smaller than those of PIC, seen in the fourth column in Table 5.So, the performance of NPIC in terms of desired mass tracking is improved with the second set of gains.The mean percentage error of the stiffness and damping for the first, second, and third SGPs obtained from both NPIC and PIC is too small, seen in the sixth and eighth columns, although the mean percentage error of position tracking is getting worse with this set of gains, seen in the second column of Table 5. Besides, the desired stiffness and damping tracking using NPIC has smaller error than those of PIC Fig. 13 Comparing the estimated impedance parameters obtained from NPIC and PIC with the desired values for the second set of gains (second SGP): the performance of NPIC is improved in terms of desired mass tracking, while the desired stiffness and damping are followed better using NPIC than PIC Fig. 14 Comparing the estimated impedance parameters obtained from NPIC and PIC with the desired values for the second set of gains (third SGP): the desired mass tracking error using NPIC is much smaller than that of PIC.The desired stiffness and damping tracking error is too small Table 5 Mean percentage error of position (in the x direction) and impedance parameters tracking: Mean (((desired value-actual value)/Desired value) *100).Gain set 1: In brief, estimating the impedance parameters and multiplying the impedance parameters error by the proportional gains K 1 , K 2 , and K 3 improve the performance of the typical PIC in terms of pose tracking and impedance model tracking and the design criteria is satisfied.Besides, the proposed controller has more freedom than the typical PIC to achieve the desired performance because of the proportional gains used in NPIC.In contrast, such gains do not exist in the typical PIC, and the performance of the PIC is not adjustable.

Conclusion
The ethical issue of using human patients for adjusting the brace, the long-time process of regulating the brace using follow-up and X-ray checks, the potential danger of repeated radiographic examinations, and passive correction pressuer control are the main disadvantages of the currently used AIS bracing treatment.Computational modeling and numerical simulations overcome all these limitations by modeling and predicting the biomechanical interaction of the AIS bracing treatment using virtual simulations before starting the brace treatment for human patients.Two challenging issues of current robotic bracing treatment of AIS are that numerical approaches, e.g., FEM have not been used for yet and analytical model of robotic brace as a single unit has not been discussed in literature.This paper proposes two computational models of the AIS bracing treatment: the MB-FE Simscape model and the analytical model.The MB Simscape model of the robotic brace is created by importing the CAD model into Simscape.An FE Simscape model of the torso is also created using the 'Reduced-order Flexible Solid' block from the 'Simscape Multibody Library'.An FE analysis is performed to compute the mechanical properties of the reduced order model.The physical interaction between the robotic brace and the torso is defined using the 'Spatial Contact Force' block from the 'Simscape Multibody Library' to form the MB-FE Simscape model.The analytical model of the AIS bracing treatment is also derived by combining the dynamics model of robotic brace as a single unit and the human torso's mass-spring-damper model.Numerical simulations and the in-vivo data from the literature are used to validate the MB-FE Simscape model and analytical model.The desired displacement of the robotic brace is applied as the primary input to both models, and the interaction force is obtained from the numerical simulations.The stiffness coefficients of the MB-FE Simscape model and analytical model are estimated using MATLAB's CFT according to the interaction force and motion measurements.The stiffness of the MB-FE Simscape model and analytical model are compared with the in-vivo stiffness data to validate both models.The in-brace correction pressure is also computed by dividing the interaction force obtained from the MB-FE Simscape model by the contact surface area.
The in-brace correction pressure is currently controlled passively by tuning the tightness of the brace's straps.An active control strategy named Novel Position-based Impedance Control (NPIC) is proposed to control the physical interaction between the robotic brace and the torso.This controller adds a feedback loop of the estimated impedance parameters to the typical PIC.A proportional controller is applied to the impedance parameters error to improve the typical PIC performance in terms of pose tracking and impedance model tracking.The significant advantage of the proposed controller is that the proportional control applied to the impedance parameters error gives more freedom to the user to tune the controller performance directly in terms of pose tracking and impedance model tracking for a given desired impedance parameters.In contrast, there is no feedback loop of the impedance parameters error in the typical PIC.The numerical MATLAB simulations are performed to verify the performance of the proposed controller in terms of pose tracking and impedance model tracking.Despite the AIS bracing treatment, the proposed controller can also be applied to other robotic applications.
In brief, two computational models were proposed to predict the biomechanical behavior of the AIS bracing treatment and a novel PIC was proposed to provide the desired interaction for AIS bracing treatment.However, further improvements are needed to improve the accuracy and performance of the proposed models and controller.One challenging issue of this study is that the interaction force obtained from the MB-FE Simscape model and the analytical model is increased (decreased) by increasing (decreasing) the stiffness and damping parameters of the 'Spatial Contact Force' block.This may increase the error between the interaction pressure reported in the literature and the interaction pressure obtained from both models.As an interesting research topic for future work, the possible solution to be proposed for this challenging issue is to create a Digital Twin (DT) of the AIS bracing treatment and use learning-based parameter identification methods and real-time force-displacement measurements to optimize the stiffness and damping parameters of the 'Spatial Contact Force' block used in the MB-FE Simscape model.The other challenging issue is that the desired mechanical impedance model for AIS treatment is unknown.The desired impedance model for AIS treatment can be determined using interaction wrench and motion measurements collected from AIS clinical treatment.Besides, the proposed controller can be powered by reinforcement learning-based algorithms to optimize the NPIC gains and overcome the negative impact of uncertainties in the dynamics model of the robotic brace and patient's torso.

Fig. 1
Fig. 1 The main contribution of the paper: a Computational modeling and predicting the physical interaction in the AIS bracing treatment, b The schematic diagram of the Novel Position-based Impedance Control (NPIC) Figure3shows the schematic diagram and the visualized form of the MB-FE Simscape model.As seen in Fig.3a, three main steps are followed to create the MB-FE Simscape model: robotic brace modeling, human torso modeling, and physical interaction definition.The first step is to create an MB Simscape model of the robotic brace using its SolidWorks model.The robotic brace was described in details in[7].The numbering of the rings, SGPs, and actuators are shown on the right-hand side of Fig.3b.To create the MB Simscape model of the robotic brace, the 'Simscape

Fig. 2
Fig. 2 Computational modeling of the AIS bracing treatment and A R D represent the rotation matrix of the frames {B}, {C} , and {D} w.r.t. the frame {A}. 1 , 2 , and 3 represent the screw angles of the frames {B}, {C} , and {D} w.r.t. the frame {A} .[s xi s yi s zi ] T , i = 1, 2, 3 denote the screw axis of the moving frames w.r.t. the fixed frame {A} .In brief, 1 = [ A P 10 T A O B T ] T , 2 = [ A P 20 T A O C T ] T , and 3 = [ A P 30 T A O D T ] T denote the pose of the first, second, and third moving rings w.r.t. the frame {A} .The pose vec- tor of the manipulator is denoted by = 1 T 2 T 3T T .The velocity vector of the manipulator χ is defined as

Fig. 4
Fig. 4 The analytical model: a The architecture of the robotic brace (the geometric parameters are only shown for the first SGP), b The schematic diagram of the mass-spring-damper model of the torso (the interaction wrench of each mass-spring-damper model at T11 level are presented by F e1 , … , F e6 ), c.The main steps to derive the ana- , and J bi3 = I 3 − A b i3× indicate the inter- mediate Jacobian matrices of the corresponding attachment points.The dynamic matrices of the moving rings are computed as follows(9)

eFig. 5
Fig. 5 The schematic diagram of the typical Position-based Impedance Control (PIC) and Novel PIC (NPIC): a The schematic diagram of the PIC, b The schematic diagram of the NPIC

Fig. 6 2 6 )
Fig. 6 Methods for numerical simulations Figure 11 represents the desired position and the position of three moving rings in the x direction obtained from NPIC and PIC.The estimated impedance parameters obtained from NPIC and PIC for three moving rings along with the desired impedance parameters are shown in Figs. 12, 13, and 14.Although the position tracking error of the NPIC is increased, the mass estimation error obtained from NPIC is too close to zero.The estimated mass obtained from NPIC is too close to the desired mass, in contrast to the estimated mass obtained from PIC.Therefore, in contrast to the PIC, the performance of the proposed controller in terms of position tracking and impedance model tracking can be improved by choosing the proper proportional gains for the NPIC.Besides, the mean percentage error of the position in the x direction and impedance parameters ( M d , C d , and K d ) obtained from NPIC and PIC for both set of gains are presented in Table 5 to quantify the comparison between NPIC and PIC.The mean percentage error is computed as follows: (((desired value − actual value)∕desired value) * 100 .The columns of Table

Fig. 11
Fig. 11 Position tracking of the moving rings in the x direction for the second set of gains (Desired-NPIC-PIC): the larger NPIC gains increase the position tracking error

Table 1
Analysis of the AIS robotic bracing treatment in terms of robot design, modeling, and control

Table 3
[13]aring the estimssated stiffness of the MB-FE Simscape model with the in-vivo stiffness data from[13]