An Overview of Mathematical Methods Applied in the Biomechanics of Foot and Ankle–Foot Orthosis Models

: Ankle–foot orthoses (AFOs) constitute medical instruments designed for patients exhibiting pathological gait patterns, notably stemming from conditions such as stroke, with the primary objective of providing support and facilitating rehabilitation. The present research endeavors to conduct a comprehensive review of extant scholarly literature focusing on mathematical techniques employed for the examination of AFO models. The overarching aim is to gain deeper insights into the biomechanical intricacies underlying these ankle–foot orthosis models from a mathematical perspective, while concurrently aiming to advance novel models within the domain. Utilizing a specified set of keywords and their configurations, a systematic search was conducted across notable academic databases, including ISI Web of Knowledge, Google Scholar, Scopus, and PubMed. Subsequently, a total of 23 articles were meticulously selected for in-depth review. These scholarly contributions collectively shed light on the utilization of nonlinear optimization techniques within the context of ankle–foot orthoses (AFOs), specifically within the framework of fully Cartesian coordinates, encompassing both kinematic and dynamic dimensions. Furthermore, an exploration of a two-degree-of-freedom AFO design tailored for robotic rehabilitation, which takes into account the interplay between foot and orthosis models, is delineated. Notably, the review article underscores the incorporation of shape memory alloy (SMA) elements in AFOs and overviews the constitutive elastic, viscoelastic, and hyperelastic models. This comprehensive synthesis of research findings stands to provide valuable insights for orthotists and engineers, enabling them to gain a mathematical understanding of the biomechanical principles underpinning AFO models and fostering the development of innovative AFO designs.


Introduction
Ankle-foot orthoses (AFOs) represent orthotic apparatus employed in the management of patients afflicted with aberrant gait patterns [1].Their primary role is to provide assistance to individuals exhibiting pathological gait and to facilitate the process of rehabilitation [2].In the context of ambulation, the term "gait cycle" is defined as the duration between two identical events commencing from the initial heel strike and culminating in the repetition of the same event for the same foot [3].The biomechanics of this movement could be a first step in the process used in the design of ankle-foot orthoses.With new technological developments, it has been observed that for professionals in engineering biomechanics, the analysis of human movement provides crucial detailed quantitative data [4].Human motion analysis involves the utilization of intricate theoretical and computational techniques aimed at elucidating the kinematic and dynamic aspects of human locomotion [5,6].Recognizing the inherent interplay between kinematics, muscle activity, and kinetic parameters, such information has been deemed indispensable for a more comprehensive comprehension of alterations in gait patterns [7].Dynamic analysis involves the comprehensive examination of motion within a multibody system, encompassing considerations of the external forces acting upon it, as well as its inherent inertial properties.A subset of these external forces is computed by solving the system's equations of motion, while others necessitate empirical acquisition [8].Within the domain of dynamic analysis, two primary problem categories emerge: forward and inverse dynamics [9].In the case of the inverse dynamic problem, it is imperative to initially gather kinematic data and ground reaction force (GRF) information [10].Forward dynamic analysis entails the simulation of the multibody system's motion over a specified time interval, integrating applied forces and initial conditions [9].In contrast, the inverse problem pertains to the determination of the dynamic forces and moments responsible for a specific motion in conjunction with the reactions occurring at each joint [11].
The resolution of inverse and forward dynamics necessitates the establishment of a biomechanical model, serving as a mathematical approximation of the actual biological system derived from external data, anthropometric parameters, and the system's equations of motion [12][13][14].A multibody system approach is employed to simulate the kinematic and dynamic environment of the system.This system comprises two or more interconnected rigid bodies with kinematic joints facilitating their motion, subjected to external forces [15].The degrees of freedom afforded by the kinematic joint configuration determine the nature of motion, introducing constraints on relative movement between the rigid bodies [16].This article's overarching objective lies in scrutinizing the existing literature to investigate mathematical methodologies applied to the analysis of ankle-foot orthoses (AFO) models.The insights gleaned from this review are poised to be of significant value to orthotists and engineers, who can leverage this information to establish a mathematical understanding of the biomechanical principles underpinning AFO models and explore avenues for the development of innovative models.

Search Strategy
The protocol type that is used in this study is the "search protocol".This protocol outlines the search strategy used to identify relevant studies for this review.It includes the databases searched, search terms used, and inclusion and exclusion criteria.
A comprehensive literature search was undertaken by querying databases such as Scopus, Google Scholar, ISI Web of Knowledge, and PubMed, spanning a time frame from 1990 to 2023.The selection of keywords employed in this systematic search process is elaborated upon below: "ankle-foot orthosis (AFO)", "numerical methods", "biomechanics", "mathematical model", "constitutive methods", "kinetic", "kinematics", and "dynamic".

Study Selection
A PRISMA (Preferred Reporting Items for Systematic reviews and Meta-Analyses) flowchart template was used for the selection of studies, and the result of this process is illustrated in Figure 1.Initially, for 196 articles, we screened titles and abstracts to identify potentially relevant articles.Subsequently, full-text reviews were conducted according to predetermined inclusion and exclusion criteria.The final selection process resulted in 42 full-text articles.
Inclusion criteria included articles in which: 1.
Studies considered for inclusion were numerical and not experimental or clinical studies alone; 2.
The means of analysis, such as equations, parameters, and methods, were clearly explained; 3.
The study was written in English and published in a peer-reviewed journal.Consistent with PRISMA guidelines, we conducted a rigorous process of data extraction.However, regarding the other PRISMA elements, such as the data synthesis and the assessment of bias, they were not relevant in our review.

Fully Cartesian Coordinates
Prior to embarking on the modeling of a system s motion within a three-dimensional space, the selection of an appropriate coordinate system becomes a pivotal prerequisite, given the array of choices available for representing a multibody system [17].This entails the utilization of unit vectors denoting joint axis rotations, presented as a column vector q, alongside a collection of points positioned at either the extremities or the articulation points within the model, collectively facilitating the representation of fully Cartesian coordinates: where V stands for vectors, P stands for points, m and n are, respectively, the number of vectors and points used to model the system.The vector composed of generalized coordinates serves as a fundamental representation of a system.The total count of these coordinates is equivalent to the number of elements within the vector, defined by the following equation: Consistent with PRISMA guidelines, we conducted a rigorous process of data extraction.However, regarding the other PRISMA elements, such as the data synthesis and the assessment of bias, they were not relevant in our review.

Fully Cartesian Coordinates
Prior to embarking on the modeling of a system's motion within a three-dimensional space, the selection of an appropriate coordinate system becomes a pivotal prerequisite, given the array of choices available for representing a multibody system [17].This entails the utilization of unit vectors denoting joint axis rotations, presented as a column vector q, alongside a collection of points positioned at either the extremities or the articulation points within the model, collectively facilitating the representation of fully Cartesian coordinates: where V stands for vectors, P stands for points, m and n are, respectively, the number of vectors and points used to model the system.
The vector composed of generalized coordinates serves as a fundamental representation of a system.The total count of these coordinates is equivalent to the number of elements within the vector, defined by the following equation: The ensemble of six coordinates comprises a trinity of orientation coordinates, which meticulously detail the spatial orientation of the system, and an additional trio of positional coordinates, meticulously pinpointing the precise location of the system concerning a static reference frame.It is of paramount importance to acknowledge that the nature of these coordinates may exhibit either a dependent or an independent character.The independent coordinates are intrinsically contingent upon the inherent degrees of freedom embedded within the system.In contrast, dependent coordinates emerge when their numerical abundance surpasses the degrees of freedom at the disposal of the system.The essence of these kinematic constraints is succinctly encapsulated within the ensuing column vector: where RBCs are the rigid body characteristics, DCs are the kinematic drivers, and JCs are the joint kinematics of the system.For a rigid body free from spatial constraints, it possesses six degrees of freedom, comprising three rotational and three translational degrees of freedom [18].The enumeration of degrees of freedom within a mechanical system corresponds to the overall count of generalized coordinates essential for elucidating the system's behavior.This tally is discerned by subtracting the number of constraints introduced by the presence of kinematic joints in the system.In essence, the degrees of freedom encapsulate the system's capacity for independent motion, while constraints restrict and define the permissible range of such motion, forming a pivotal aspect of mechanical system analysis [19].

Kinematic Analysis
This analysis pertains to the study of motion, irrespective of the forces instigating it, and is sometimes referred to as the "initial position problem".The nomenclature aligns with the task of determining the positions of all components within the system, factoring in the input provided by driving constraints [20].From a mathematical perspective, this entails determining the generalized coordinates (q) that satisfy the set of kinematic constraint equations at each time step of the analysis, described as follows: To solve Equation (3), the Newton-Raphson method (NRN) is employed, known for its quadratic convergence properties near the solution [21], achieved by linearizing the system at q i : Φ(q, t) ∼ = Φ(q i ) + Φ q (q i )(q The resolution for the partial derivatives of each kinematic constraint concerning the vector of generalized coordinates finds its representation in the form of the Jacobian matrix of constraints, denoted as Φ q (q i ), as stipulated in Equation (6): where nc denotes the number of dependent coordinates and nh signifies the count of constrained equations.The vector representing an approximation of the correct system position is derived as the solution to Equation ( 4).An iterative equation, denoted as q i+1 , is established as follows: Φ(q i ) + Φ q (q i )(q i+1 − q i ) = 0.
To determine the system's velocity, one simply needs to differentiate the kinematic constraint Equation (4) with respect to time, which results in the following expression: where ∂Φ(q,t) ∂t represents the vector v(t), consisting of the partial derivatives of Φ with respect to time; ∂Φ(q,t) ∂q signifies the Jacobian matrix, and ∂q ∂t represents the vector of generalized velocities, .q. Thus, the velocity equation is articulated as: Differentiating the constraint velocities from Equation ( 8) with respect to time yields the vector of constraint accelerations: .. Φ(q, .q, ..
Since the vector γ q, .q, t is defined as: then, the acceleration equation results in:

Dynamic Analysis
To tackle the dynamic analysis problem for a multibody system, it is imperative to possess a comprehensive understanding of the external inertia and forces acting upon each individual rigid body within the system [22].Several methods can be employed to derive the equations of motion, with one of them being the principle of virtual power, as presented below.In this principle, the sum of the virtual power generated by the external forces at each moment equates to zero [16]: where, the vector f, representing all the forces generating virtual power, is expressed as: In the context of the equation set, M ..
q represents the vector signifying the inertial forces, where M denotes the comprehensive global mass matrix and .. q symbolizes the vector delineating generalized accelerations.Simultaneously, g stands as the vector encompassing generalized applied forces.It is noteworthy that the internal forces, organized in actionreaction pairs, do not contribute to the concept of virtual power and, therefore, are not factored into Equation ( 13).Nevertheless, their determination is attainable through the utilization of the Lagrange multiplier technique, presented as follows: In this equation, g Φ designates the vector characterizing generalized forces, while λ constitutes the column vector housing the Lagrange multipliers.The magnitude of these internal constraint forces is ascertained through these multipliers, whereas their specific direction is contingent on the information encoded in the Jacobian matrix.The equation encapsulating the concept of virtual power is deduced by amalgamating the expressions found in Equations ( 13)-( 15): The holistic equation of motion, presented comprehensively below, serves as the cornerstone to unearth the system's unknown variables: This method proves to be instrumental in bringing together and solving the equations of motion ( 17) for the biomechanical system.Depending on the selected approach for solving these equations, researchers can engage in either forward dynamic analysis (Simulation) or inverse dynamic analysis.The latter, in particular, plays a pivotal role in deducing the external forces governing the observed motion, thereby unveiling the hitherto unknown external and internal forces.This analytical approach yields invaluable results, significantly contributing to the field [23].

Two-Degree-of-Freedom AFO for Robotic Rehabilitation
The capability of accommodating eversion-inversion motion in ankle-foot orthoses is contingent on the flexibility of the material employed.However, constraints on the normal eversion-inversion range contribute to discomfort and hinder the emulation of a natural ankle motion [24].In this context, the discussion centers on an ankle-foot orthosis (AFO) endowed with two degrees of freedom.These degrees of freedom encompass plantarflexiondorsiflexion and eversion-inversion motions.Notably, the plantarflexion-dorsiflexion motion is actively regulated by a DC servomotor, whereas the eversion-inversion joint operates passively, facilitated by a damper and a torsion spring.The latter components are instrumental in restricting the range of motion beyond a specific threshold by applying spring force [25].

Model of the Foot
The intricate motion of the foot is a multifaceted process that transpires across three distinct planes and revolves around three separate axes [26].The orientations of the remaining two joint axes and the mapping of these joint axes onto the inertial frame are integral components of this intricate movement and are depicted in Figure 2. The inertial frame is affixed to the axis and the shank, ensuring a stable reference frame [27].
It is pertinent to note that → Z 0 symbolizes the knee joint axis, → Z 1 denotes the plantarflexiondorsiflexion joint axis, and → Z 2 represents the eversion-inversion joint axis.These axes are considered when the foot is in firm contact with the ground and the shank is perpendicular to the ground, thus setting the reference projection.Sequential rotations are essential to transition

Model of the Foot
The intricate motion of the foot is a multifaceted process that transpires across three distinct planes and revolves around three separate axes [26].The orientations of the remaining two joint axes and the mapping of these joint axes onto the inertial frame are integral components of this intricate movement and are depicted in Figure 2. The inertial frame is affixed to the axis and the shank, ensuring a stable reference frame [27].It is pertinent to note that  ⃗ symbolizes the knee joint axis,  ⃗ denotes the plantarflexion-dorsiflexion joint axis, and  ⃗ represents the eversion-inversion joint axis.These axes are considered when the foot is in firm contact with the ground and the shank is perpendicular to the ground, thus setting the reference projection.Sequential rotations are essential to transition from  ⃗ to  ⃗ and subsequently from  ⃗ to  ⃗ , facilitating the determination of rotation matrices required to reach the desired frame.The representation relies on spherical coordinates, which are inherently linked to Cartesian coordinates as follows: [x, y, z] = [rcos γsin ϕ, rsin γsin ϕ, rcos ϕ]. ( The three Z vectors within the inertial frame can be expressed as follows: The angles γ and ϕ can be calculated based on the projection of joint axes in the inertial frame.The projections of the → Z i axis in the inertial frame are represented as (X i ′ , Y i ′ , Z i ′ ) for i = 1, 2, and can be determined using the following equations: Once the orientations of The Denavit-Hartenberg parameters θ i , α i , a i , and d i are then calculated for the transformation of frames f i to frame f i+1 [28].The transformation matrix can be represented as: The position loop closure equation is used to find the a i and d i : In this equation, T 3 0 is computed independently using The parameters a 0 , d 0 , d 1 , a 2 , and d 2 are assumed based on the foot's geometry, while a 1 , a 3 , and d 3 are calculated from three equations obtained from the previous loop closure equation.The complete kinematic model becomes known once all these parameters are determined.
Leveraging the assistance of the kinematic model, one can proceed to formulate the dynamic model, which, in turn, streamlines the process of designing the orthosis.The dynamic model of the system operates under the framework of the Lagrangian formulation and assumes a fixed shank, simplifying the overall analysis.This formulation is derived from the potential and kinetic energies of the system.As additional springs and masses are integrated into the system, the dynamic model must be adjusted to accurately represent the potential and kinetic energies of the system.The equations of motion are expressed as: Here, D(q) denotes the (2 × 2) inertia matrix, C q, .q represents the (2 × 2) terms for Coriolis and centrifugal effects, G(q) stands for the (2 × 1) vector of gravity torque, and q = [θ 1 , θ 2 ] T is the vector encompassing all joint variables.Once the complete motion of the system is known through the Newton-Euler equations, and the force-torque sensors provide

Model of the Orthosis
, the remaining unknowns can be uniquely determined.The set-point control of the ankle-foot orthosis is a noteworthy aspect of this underactuated system, characterized by the presence of a damper and a spring at the passive eversion-inversion joint, is represented in Equation (28) and can be expressed as: thy aspect of this underactuated system, characterized by the presence of a damper and a spring at the passive eversion-inversion joint, is represented in Equation ( 28) and can be expressed as: Therefore, it is imperative to verify the stability of the control law for the system.This is achieved by defining a Lyapunov function, VL, as the total energy of the system, along with an additional error term: Within this equation,  represents the potential energy attributed to gravity,  is a positive constant, and  =  −  signifies the error associated with the first joint.Moreover,  is another positive constant ensuring that  > 0 remains greater than zero.The primary aim here is to establish that throughout any trajectory of the system, the function  steadily decreases, eventually converging to zero.This phenomenon implies that the machine is progressively advancing towards the desired configuration denoted Therefore, it is imperative to verify the stability of the control law for the system.This is achieved by defining a Lyapunov function, V L , as the total energy of the system, along with an additional error term: Within this equation, V m represents the potential energy attributed to gravity, K p is a positive constant, and e 1 = θ 1 − θ 1d signifies the error associated with the first joint.Moreover, V 0 is another positive constant ensuring that V L > 0 remains greater than zero.The primary aim here is to establish that throughout any trajectory of the system, the function V L steadily decreases, eventually converging to zero.This phenomenon implies that the machine is progressively advancing towards the desired configuration denoted as q d .Assuming a constant q d and making use of the properties of a skew-symmetric matrix, specifically, .q T .D − 2C .q, being equal to zero, the application of Equations ( 29) and ( 30), in conjunction with the time derivative of V L , leads to the formulation of the ensuing equation: .
Here, K td is the damping coefficient.By selecting a PD control law = −K p e 1 − K d .θ 1 : .
The analysis demonstrates that V L decreases as long as K p can be selected to be sufficiently large to drive e 1 close to zero, causing .θ 1 to approach the desired set point.
The control torque is provided by the actuator within the simulation, although it is important to acknowledge that, in practical human use, an internal torque will also be generated by the muscles.In such instances, the net torque required by the actuator is determined by the equation: Here, τ M represents the machine component of the torque, while τ M signifies the human component.τ M can be ascertained in real-time using force-torque sensor data.Once τ M is known, it allows for the calculation of the necessary motor torque based on Equation (34).

SMA-Element-Based AFO
Various methodologies have been introduced to model the behavior of shape memory alloys (SMAs).One of the models, proposed by Brinson and Huang, segregates the martensite volume (ξ) into two distinct components: temperature-induced martensitic (ξ T ) and stress-induced martensitic (ξ S ) [29].To delineate ξ S and ξ T in the context of the shape memory effect (SME) at low-temperature M s , the Liang and Roger transformation equation is employed [30].The division of the martensite volume fraction (ξ) is expressed as follows [31]: The constitutive equation that relates stress (σ), strain (ε), and temperature (T) in terms of the martensitic volume fraction (ξ) is given by: Here, ε L denotes the maximum strain that can be recovered [31], while E represents the elastic modulus [30]: where E M stands for the elastic modulus in the martensite phase, and E A corresponds to the austenite phase.The phase transformation is initiated at specific temperatures, known as M s (martensite start), A s (austenite start), M f (martensite finish), and A f (austenite finish) temperatures under zero-stress conditions.The stresses associated with the detwinning start and finish processes are symbolized as σ cr s and σ cr f , respectively.C M signifies the constant that illustrates the influence of stress coefficients, and C A (as depicted in Figure 4) represents the impact of stress on the transformation temperatures of the martensite and austenite phases [30].To determine the martensitic fractions concerning both temperature and stress (as depicted in Figure 4), the progression of equations can be delineated as follows: For   and  +  ( −  )   +  ( −  ), To determine the martensitic fractions concerning both temperature and stress (as depicted in Figure 4), the progression of equations can be delineated as follows: For T > M S and Also, for T < M S and σ cr s < σ < σ cr f , where, if M f < T < M S and T < T 0 , Here, the parameters (T 0 , ξ S 0 , ξ T 0 , and ξ 0 ) denote the initial state of the material.Additionally, a A and a M are material constants defined as: The behavior of the shape memory alloy is characterized by a nonlinear function involving temperature, stress, and pressure, which are considered to be independent and unrelated to each other [32].
Moreover, concerning the mathematical techniques and models within the framework of SMA (shape memory alloy), the fundamental framework of the model was elucidated through the utilization of a rheological model [33].The foundational kinematic assumptions of the finite deformation model and the Helmholtz free energy function were detailed [34], and the exploitation of the Clausius-Duhem inequality was employed to derive the constitutive inequalities [33].These inequalities serve as the basis for establishing the evolution equations.

Linear Elastic Constitutive Model
The linear elastic constitutive model, deeply rooted in Hooke's Law, a fundamental principle in mechanics, has had a profound impact on various fields of science and engineering.This model establishes a proportional relationship between the strain (change in shape or size) in a material and the applied stress (force) within the elastic limit of that material [35].This linear relationship forms the cornerstone of understanding material behavior, and it has proven to be invaluable in the context of ankle-foot orthosis (AFO) design.To understand the significance of the linear elastic constitutive model, we must first trace its origins to Hooke's Law, a cornerstone principle in mechanics [36].In mathematical terms, this relationship is expressed as: In this equation, σ represents the stress applied to the material, E stands for the material's Young's modulus (a constant that characterizes the material's stiffness), and ϵ denotes the strain.
This fundamental law is a precursor to the linear elastic constitutive model and serves as its conceptual foundation [37].Hooke's Law established the initial understanding that the response of a material to an applied force is linear within the elastic limit, and this notion is encapsulated in the linear elastic constitutive model.
The linear elastic constitutive model is often presented in terms of tensorial quantities, particularly emphasizing the fourth-order elasticity or stiffness tensor.This tensor is a comprehensive mathematical framework that describes the relationship between stress and strain in a material.The generalized Hooke's Law, which encapsulates this relationship, defines a linear relation among all the components of the stress and strain tensor and is represented as [35]: In this equation, σ ij represents the components of the stress tensor, ϵ kl denotes the components of the strain tensor, and C ijkl signifies the components of the fourth-order stiffness tensor.For a three-dimensional context, this tensor comprises a total of 81 components, while two-dimensional problems reduce the count to 16 components.The inclusion of tensorial quantities in this model underscores its versatility and applicability to a wide range of materials and structures.This tensor-based approach is capable of describing complex mechanical behaviors and allows for the modeling of materials with varying properties in different directions.This mathematical framework is particularly well-suited for materials with anisotropic properties, where stiffness and elasticity differ along different axes [38,39].
The utility of the linear elastic constitutive model extends far beyond theoretical considerations.This model has found widespread use in the field of biomechanics, where understanding the mechanical properties of human tissues, as well as the design and analysis of orthotic devices like AFOs is of paramount importance.Researchers have harnessed this model to gain deeper insights into the behavior of human tissues, particularly in the context of AFO design.For example, studies have employed finite element analysis, a computational technique for solving complex problems, to investigate the three-dimensional behavior of the human foot during activities like standing.These studies have employed the linear elastic constitutive model to describe the mechanical properties of the foot tissues, offering valuable insights into the stresses and strains experienced by the foot under various conditions [37].
Similarly, computational models have been developed to study the effect of foot constraint on ankle injuries.These models incorporate the linear elastic constitutive model to represent the mechanical behavior of the ligaments and foot structures, providing a deeper understanding of how these structures respond to various loads and constraints [40].Moreover, the linear elastic constitutive model has been instrumental in the design and optimization of AFOs.Researchers and engineers have utilized this model in virtual prototyping processes for the manufacture of passive-dynamic AFOs [41].By simulating the bending stiffness of the orthotic device, they can optimize its design for improved performance and patient comfort.This process has accelerated the development and refinement of AFOs, leading to more effective orthotic devices [42].
One of the most significant contributions of the linear elastic constitutive model to AFO design is the ability to customize orthotic devices to cater to the unique biomechanical needs of individual patients.This is especially critical in the field of orthotics and rehabilitation, where every patient's condition and requirements are distinct [43].The model's versatility is apparent in its influence on material selection for AFOs.The choice of materials is a fundamental aspect of orthotic design, as it directly impacts the device's performance.The linear elastic constitutive model offers orthotists a starting point for material selection.Stiffer materials, as guided by the model, are preferred when the goal is to enhance support and stability in AFOs.These materials are well-suited for patients requiring significant structural support or conditions requiring rigid bracing [44].
Conversely, more flexible materials are advantageous when the objective is controlled deformation, allowing for the mimicry of natural foot movement [45].These materials are particularly useful in cases where patients need a more dynamic and accommodating orthotic device.The model empowers orthotists to tailor AFOs to the specific biomechanical needs of individual patients, ensuring an optimal balance between support and comfort.This customization is vital for enhancing patient compliance and overall quality of life [46].The ability to select and customize materials based on the linear elastic constitutive model offers orthotists a systematic and science-based approach to AFO design.This approach ensures that the materials used are not only appropriate for the intended purpose but also contribute to the overall effectiveness of the orthotic device.
The linear elastic constitutive model offers several distinct advantages that have made it a vital tool in the field of AFO design and biomechanics [47]: (1) simplicity-The model's linear approach is straightforward and relatively easy to understand, providing a clear foundation for analyzing material behavior within the elastic limit; (2) predictability-the model's linear nature allows for predictable and repeatable results, making it a reliable tool for engineers and orthotists; (3) versatility-the inclusion of tensorial quantities in the model enables it to describe a wide range of materials, including anisotropic materials with varying properties along different axes; (4) customization-the model empowers orthotists to tailor AFOs to individual patient's needs, balancing support and comfort effectively; (5) accelerated development-the use of this model in virtual prototyping processes has accelerated the development and refinement of AFOs, leading to more effective orthotic devices.
However, like any model, the linear elastic constitutive model has its limitations [48,49]: (1) linear elastic assumption-the model simplifies material behavior to linear elasticity.In reality, many materials exhibit nonlinear behavior under certain conditions, such as plastic deformation or material failure.The model cannot accurately represent such phenomena; (2) real-world complexities-designers and orthotists must use this model judiciously, considering the specific properties of materials and the real-world conditions to which AFOs will be subjected.In practice, materials often exhibit complex behaviors beyond the scope of linear elasticity; (3) limited to elastic range-the model is applicable only within the elastic limit of materials.Beyond this limit, materials may exhibit irreversible deformations and complex behavior that cannot be accurately predicted by the linear model.
In practice, it is essential to acknowledge these limitations and apply the model within its appropriate scope.Orthotists and engineers should use the model as a foundational tool but be prepared to account for nonlinear behaviors and other complexities when designing and manufacturing AFOs for real-world applications [50].

Viscoelastic Constitutive Model
Viscoelasticity is a fascinating property observed in materials that exhibit both viscous and elastic characteristics when subjected to deformation.This unique behavior, where materials display both time-dependent strain and stress responses, has significant relevance in various fields [51], including material science, engineering, and particularly in the construction of ankle-foot orthoses (AFOs) [52].The viscoelastic constitutive model, a mathematical model designed to characterize and predict the mechanical behavior of viscoelastic materials, plays a pivotal role in understanding and utilizing these materials effectively [53].
The mechanical behavior of viscoelastic materials is inherently complex due to its timedependent nature [54].To accurately capture this behavior, the viscoelastic constitutive model employs intricate mathematical representations.One common form of this equation can be expressed as: In this equation, σ(t) represents the stress at time t, ϵ(τ) denotes the strain at time τ, and G(t − τ) is a relaxation function that describes how the material response changes over time.This equation captures the complex interplay between stress and strain in viscoelastic materials and has a direct application in the design of AFOs.
Viscoelastic materials exhibit several distinctive characteristics that set them apart from purely elastic materials.Understanding these properties is fundamental to comprehending the behavior of these materials in real-world applications: (1) hysteresis-viscoelastic materials exhibit hysteresis, a fundamental property characterized by the energy dissipation within the material during loading and unloading cycles.This behavior underscores the role of viscous properties in dissipating energy, a critical aspect of viscoelasticity [55]; (2) creep-creep is another defining characteristic of viscoelastic materials.It refers to the increase in strain under constant stress over time.This phenomenon is often observed in materials used in AFOs, reflecting the materials' tendency to deform under prolonged stress, making it particularly relevant in orthotic design; (3) stress relaxation-stress relaxation represents the decrease in stress under constant strain over time, illustrating the material's ability to adapt its response to changing conditions.This property is crucial in AFO design, as orthotic devices must provide support while accommodating changes in stress over time, especially during activities like walking [56].
The viscoelastic constitutive model, while offering a more accurate description of material behavior, introduces increased complexity to the analysis.The time-dependent nature of the material response necessitates the use of more advanced numerical methods for solving problems involving viscoelastic materials.This complexity can involve intricate differential equations, numerical integration techniques, and a deeper understanding of the underlying material properties [57].
The unique properties of viscoelastic materials have a direct application in the design and optimization of AFOs, which must accommodate the dynamic and varied loads experienced during activities such as gait, all while ensuring patient comfort and support.By applying the viscoelastic constitutive model, orthotists and engineers can predict how viscoelastic materials will deform and respond to diverse loading conditions.This approach allows for the creation of more effective AFOs tailored to individual patient needs, optimizing the balance between support and comfort [58].
Furthermore, the understanding of viscoelastic material behavior enables AFO designers to consider time-dependent effects in orthotic design, enhancing the overall performance of these devices.It ensures that AFOs are not solely designed based on the assumption of linear elasticity, which is a simplification often used for purely elastic materials.By accounting for viscoelasticity, orthotic devices can more accurately mimic natural biomechanics, providing patients with a higher level of comfort and functionality during daily activities.

Hyperelastic Constitutive Model
Hyperelastic constitutive models, often referred to as strain energy density functions or material models, are fundamental tools in the field of biomechanics [59] and material science [60].These models, rooted in the theory of elasticity, provide a systematic approach to characterize the behavior of materials, particularly those used in ankle-foot orthoses (AFO).Hyperelastic models are especially valued for their ability to capture the nonlinear, timeindependent response of materials to deformation.One of the most recognized hyperelastic models is the neo-Hookean model, exemplifying the essence of these equations [53].
At the heart of the hyperelastic constitutive model lies the mathematical representation that characterizes the material behavior [61].This representation hinges on the definition of a strain energy density function, W, in terms of the deformation gradient, F. The neo-Hookean model, an archetype of hyperelastic models, can be expressed as [62]: Here, W symbolizes the strain energy density function, while µ and λ are material constants.I 1 represents the first invariant of the deformation tensor, and J denotes the determinant of the deformation gradient.The neo-Hookean model captures the essential characteristics of hyperelastic materials, providing a foundation for understanding their complex behavior [63].
The hyperelastic constitutive model plays a pivotal role in AFO design.AFOs must accommodate the dynamic and varied loads experienced during gait while ensuring patient comfort and support [64].The application of hyperelastic models enables orthotists and engineers to predict the deformation of materials under diverse loading conditions, thereby ensuring AFOs deliver the requisite support and flexibility to patients with a wide spectrum of pathologies [65].
The hyperelastic constitutive model offers several advantages in the realm of AFO design [66], including: (1) nonlinearity-these models capture the nonlinear behavior of materials, providing a more accurate representation of the intricate deformations occurring in AFOs; (2) customization-orthotists can tailor AFO designs by adjusting material parameters, and matching the orthotic device to the specific biomechanical needs of individual patients; (3) predictive modeling-the models allow orthotists and designers to predict material responses under various loading conditions, facilitating the creation of more effective and patient-specific AFOs.Despite their advantages, the hyperelastic constitutive model is not without limitations [67].These include: (1) material complexity-not all materials can be accurately described by a single hyperelastic model.More complex materials may require the incorporation of multi-physics models to account for their unique characteristics; (2) parameter estimation-accurately determining material parameters for hyperelastic models can be a challenging endeavor, often requiring extensive experimentation and testing.
The hyperelastic constitutive model, with the neo-Hookean model as a prominent example, plays a crucial role in AFO design [68].Their ability to accurately capture nonlinear material behavior, enable customization, and support predictive modeling makes them indispensable tools in the creation of orthotic devices that enhance patient mobility, comfort, and overall quality of life.The broader implications of these models extend into material science, biomechanics, and personalized healthcare, where the principles of hyperelasticity continue to shape how we understand and manipulate the physical world [69].

Conclusions
In summary, this review of the literature on mathematical methods applications in the biomechanics of foot and ankle-foot orthosis models has uncovered several crucial insights, shaping the future of this field.Regarding the nonlinear optimization and Cartesian coordinates, it is found that employing nonlinear optimization in fully Cartesian coordinates enhances our understanding of AFOs' kinematic and dynamic aspects.This opens avenues for more advanced and effective AFO designs.In terms of robotic rehabilitation with two-degree-of-freedom AFO, the exploration of a two-degree-of-freedom AFO for robotic rehabilitation stands out as a promising approach.Considering the interplay between foot and orthosis models, this innovation has the potential to significantly impact rehabilitation strategies.Regarding the shape memory alloy (SMA) elements in AFOs, the review underscores the growing trend of incorporating shape memory alloy (SMA) elements in AFOs.This suggests exciting possibilities for further research and development, offering potential benefits in AFO functionality.Additionally, the diversity of constitutive models, including linear elastic, viscoelastic, and hyperelastic models, emerged as a key finding.Understanding these models is crucial for tailoring AFO designs to individual patient needs, providing a more personalized approach to orthotic interventions.
Building on these insights, future research in AFOs should prioritize: (1) innovative mathematical models-explore new mathematical models and optimization techniques to enhance precision in AFO design; (2) clinical implications of novel designs-investigate the long-term clinical implications and patient outcomes associated with innovative AFO designs, especially those incorporating SMA elements; (3) interdisciplinary collaborationencourage collaboration between orthotists, engineers, and clinicians to bridge the gap between theoretical advancements and practical applications in patient care.

Figure 1 .
Figure 1.Selection of studies based on the PRISMA flow diagram.

Figure 1 .
Figure 1.Selection of studies based on the PRISMA flow diagram.

Z 2 ,
facilitating the determination of rotation matrices required to reach the desired frame.

Figure 2 .
Figure 2. Illustration of the orientations of the plantarflexion-dorsiflexion joint axis ( ⃗ ) and the eversion-inversion joint axis ( ⃗ ), along with their corresponding kinematic representations [25], with permission from the American Society of Mechanical Engineers ASME (License Number: 1408140-1).

Figure 2 .
Figure 2. Illustration of the orientations of the plantarflexion-dorsiflexion joint axis (

→ Z 1
and → Z 2 in the inertial frame are established, the Denavit-Hartenberg parameters can be computed, allowing for the formulation of rotation matrices by defining the →

Figure 3
Figure3depicts the free-body diagram of all the components within the system.The machine segments are represented by rectangular shapes, while the human segments are denoted by bony and oval shapes.The weight of each body acts at its center of mass.Sensor forces and torques are illustrated by circles.The subscripts M correspond to the machine components, H to the human components, and S signifies the sensor component.

Figure 3 .
Figure 3. Machine and human parts of the shank and ankle [25], with permission from the American Society of Mechanical Engineers ASME (License Number: 1408141-1).

Figure 3 .
Figure 3. Machine and human parts of the shank and ankle [25], with permission from the American Society of Mechanical Engineers ASME (License Number: 1408141-1).