Qualitative Analysis of Solutions of a Nonlinear Biomechanical Model

Ifeoma D. Omoko^1*, Ayinla A. Abdurasid², Jamiu R. Akewushola¹, Ola A. Salami³

¹Department of Mathematics, Lagos State University, Lagos, Nigeria.
²Department of Mathematics, Lagos State University of Science and Technology, Lagos, Nigeria.
³Department of Mathematics, Federal College of Fisheries and Maritime Technology, Lagos, Nigeria.

DOI: 10.4236/jamp.2024.1211242                

Abstract

The behavior of objects in motion is described by the equations of motion, which are basic concepts in mathematical physics. These equations are useful in explaining how forces and torques cause body components to move around a joint when applied to joint movement, especially in biomechanics. In the orthopedic industries, biomechanics is widely used to develop orthopaedic implants for a range of body joints, dental parts, external fixations, exoskeletons, and other medical uses. In this case, the motion of a phenomenon is described using a nonlinear differential model. One of the most effective approaches for describing the qualitative behavior of a dynamical system is the introduction of Lyapunov methods. Stability analysis and boundedness of solutions of a nonlinear differential equation model, particularly in the context of knee joint movement, entails analyzing how minor perturbations (like changes in force, position, or velocity) influence the behavior of the joint and remain within a finite range over time, respectively. The goal is to determine whether the system returns to a steady state (stable) or becomes unstable when subjected to these small changes. The effect of viscous damping, external input, and angular motions at different times in seconds are all controlled to govern the shank knee movement surrounding the knee joint. Numerical simulations with Matlab and Mathematica are drawn to demonstrate the effectiveness of the shank knee motion around the knee joint.

Keywords

Nonlinear Differential Equations, Biomechanics, Stability, Boundedness, Lyapunov Methods, Knee Joint

1. Introduction

Engaging in regular physical activity is one of the most significant things people can do to decompress and relieve themselves of all mental tension. Overuse, direct impact, or applying force beyond what a bodily part is anatomically capable of withstanding are the main causes of sports injuries, according to Better Health Channel ([1] [2]). As the largest and most intricate joint that is an essential component of the human body, the knee joint is also the one that is most vulnerable to damage. Research on knee biomechanics is concerned with knee mobility, joint contact, soft tissue deformation, and the impact of force on joints. Because knee movement, joint contact, and soft tissue deformation are the main topics of study in knee biomechanics, it is crucial to understand how force affects joint motion. Predicting stresses in various knee segments during daily activities is much easier with an understanding of biomechanics. The maximum compressive stress of the fully extended knee during the gait cycle, contact pressure, contact area, and other biomechanical parameters were assessed, Wang ([3] [4]).

Biomechanical analysis shows from the results of some researchers ([5]-[7]) that there are two major courses, namely Total Knee Replacement (TKR) and Active Orthosis. Total Knee Replacement (TKR) is unavoidably conducted for many patients with severe knee injuries. According to the Instanbul Medical Association, Türkiye [4], “When knee pain becomes a barrier to your daily life, knee replacement surgery emerges as the best solution. This procedure aims to remove pain and reinstate the normal functions of your knees. Surgeons skillfully remove damaged bone and cartilage from the thighbone, shinbone, and kneecap, paving the way for a renewed lease on mobility. Then, they replace it with an artificial joint (prosthesis) crafted from metal, high-grade plastics, and polymers. Surgeons assess your knee’s motion, stability, and strength to determine if a knee replacement is necessary. Additionally, they utilize X-rays to evaluate the extent of the damage”.

Surgeons skillfully consider a range of prostheses and surgical techniques, such as age, weight, lifestyle, knee size and shape, and overall health, which are crucial in determining the most suitable approach. During knee replacement surgery, damaged knee components are substituted with metal and plastic counterparts [2].

An orthosis, by the British Association of Prosthetists and Orthotists, is an externally applied device used to modify the structural and functional characteristics of the neuro-muscular and skeletal systems. Examples of orthosis include insoles; footwear and adaptations; spinal braces, callipers; knee, elbow, shoulder, wrist and hand braces [1]. They produce the joint torques needed both to compensate for gravity and loading forces and to recreate natural human movements.

Research on the nonlinear model describing the motion of the shank around the knee joint is quite numerous. For example, Amélie Chevalier et al. [6] described an internal model control for shank movement around the knee joint as

$L_T\frac{{\text{d}}^2\theta \left(t\right)}{\text{d}{t}^2}+B\frac{\text{d}\theta \left(t\right)}{\text{d}t}+K\theta \left(t\right)+M_{T}g{r}_T\sin \left(\theta \left(t\right)\right)=T\left(t\right)$ (1)

where $L_T$ length of the shank, $r_T$ which is the distance between the center of mass of the shank and the knee joint axis, $m_T$ which is the mass of the shank, g which is the gravitational constant.

Amélie Chevalier et al. [6] extend their results [5] to nonlinear predictive control of shank movement around the knee joint using the equation of motion (1),

Other findings on Biomechanical analysis include Maquet [8], Winter [9], Zavatsky [10].

2. Statement of the Problem

In the study of differential equations, control systems, etc, stability problems and boundedness can be addressed using the Lyapunov Methods. The stability theory has produced amazing discoveries and has a broad application in concrete problems of the real world. Ultimately, a system’s ability to maintain stability is one of its most crucial requirements, and without it, the system would be ineffective and possibly chaotic [11].

From the results of Amélie Chevalier et al. ([5] [6]), Sukhman Kaur [7], replacing $K\theta \left(t\right)+M_{T}g{r}_T\sin \left(\theta \left(t\right)\right)$ of Equation (1) by a continuous functions, $bG\left(x\right)\left(t\right)$ , this mechanical model is described by nonlinear differential equations,

$\ddot{x}+a\dot{x}+bG\left(x\right)+c\sin \left(x\right)=\delta T\left(t,x,\dot{x}\right)$ (2)

where $G\left(x\right)\left(t\right)$ is a continuous function.

This study presents a qualitative analysis in terms of stability and boundedness of solutions of a nonlinear biomechanical model describing the movement of the shank, around the knee joint by the Lyapunov Methods. Other results on Lyapunov Methods include Olutimo and Omoko ([12]-[15]).

3. Mathematical Modelling

A nonlinear differential equation of motion model can be highly effective in addressing knee joint injuries by offering a precise and in-depth depiction of the joint’s behavior under various conditions. This model captures a more accurate view of these dynamics, capturing the essential non linearities in stiffness, damping, and muscle force generation. The knee joint covers elaborate interactions between bones, cartilage, ligaments, tendons, and muscles. This model can represent these interactions more accurately than linear models. In fact, this model can be used to design more effective knee braces and supports. These devices can be tailored to provide the right amount of support and limit harmful motions without unduly limiting normal movement. The model can simulate how different brace designs affect knee mechanics under various conditions, ensuring that the brace helps prevent injury without compromising performance. From the system

$\ddot{x}+a\dot{x}+bG\left(x\right)+c\sin \left(x\right)=\delta T\left(t,x,\dot{x}\right)$

$G\left(x\right)\left(t\right)$ , a continuous function, the parameters are oftentimes expressed as;

$\ddot{x}=\frac{\text{d}{x}^2}{\text{d}{t}^2}$, the angular acceleration;

$\dot{x}=\frac{\text{d}x}{\text{d}t}$, the angular velocity, x = output angle;

$a=\frac{B}{I}$ , ratio of viscous damping to the inertia;

$b=\frac{K}{I}$ , ratio of the knee joint stiffness to the inertia, $c=Mg$ .

Unlike the previously reported results on Shank knee models, a complete Lyapunov function will be introduced to analyze the stability and boundedness of the output angles.

4. Stability Analysis

Stability analysis determines whether, after a small perturbation, the knee joint will return to its equilibrium position (stable), oscillate around it (marginally stable), or move away from it (unstable). Using a nonlinear equation of motion, we are to apply the Lyapunov methods to assess the motion of the Shank knee around the knee joint. A scalar function, say $V\left(X\left(t\right)\right)$ , will be presented as the Lyapunov function, with its positive definiteness and the requirement for a negative definite derivative, $\dot{V}\left(X\left(t\right)\right)$ along the trajectory of the system examined. The consideration of stability property on the angles in the equation of motion occurs when these conditions are met. Subsequently, we propose some basic concepts and a theorem to support the stability analysis of the shank knee motion in the following section.

The zero solution of Equation (2) is STABLE if given $ϵ>0$ , $t_0>0$ , there exists a $\delta \left(t_0,ϵ\right)>0$ such that whenever,

$\left|X_0\right|<\delta \left(t_0,ϵ\right)$ , $\left|X\left(t,t_0\right)\right|<ϵ$ for all $t\ge t_0$ .

The zero solution of Equation (2) is ASYMPTOTICALLY STABLE (AS), if it is stable and in addition, there exists $\alpha \in \left[t_1,t_2\right]$ , $t_0\le t_1\le t_2\le t$ such that if $X_0\le \delta \left(t_0,\alpha \right)$ , then

$\left|X\left(t;t_0,X_0\right)\right|\Rightarrow 0$ as $T\to \infty$ .

Drawing from these two definitions, we present the assumptions that stability theory can be proved.

Assumptions

In addition to the assumptions imposed on the model parameters, a, b, c, $L_T,g,K,m,\cdots$ in (1) as non-negative values, this study also assumed that

• $\left|\Psi \left(x\right)\right|>b$

• $ab>0$

To establish our result, Equation (2) is linearized in the absence of $\delta T\left(t,x\right)$ to yield the form

${\dot{x}}_1=x_2$

${\dot{x}}_2=a\dot{x}-bG\left(x\right)-c\sin \left(x\right)$

and introducing $x_1\Psi \left(x_1\right)=bG\left(x_1\right)+c\sin \left(x_1\right)$ , we select a Lyapunov function $V\left(x\left(t\right)\right)$ defined by:

$2V\left(x_1,x_2\right)={\left(a{x}_1^2+x_2\right)}^2+x_2^2+4{\displaystyle {\int }_0^{x_1}\Psi \left(s\right)\text{d}s}$ (3)

From Equation (3), it is obvious that $V\left(0,0\right)=0$ , and V is positive definite. That is,

$V\left(x_1,x_2\right)\ge x_2^2+{\left(a{x}_1+x_2\right)}^2+x_1^2\Psi \left(x_1\right)$

$V\left(x_1,x_2\right)\ge x_2^2+{\left(a{x}_1+x_2\right)}^2+b{x}_1^2$

Equation (2), therefore, satisfies the condition

$V\left(x_1,x_2\right)\ge u_1\left(x_1^2+x_2^2\right)$ (4)

where $u_1$ is a positive constant, provided $\left|\Psi \left(x_1\right)\right|>b$ , $ab>0$ .

The time derivative of $V\left(x_1,x_2\right)$ along the trajectory (2) is thus simplified using the notation as:

$\dot{V}=\frac{\partial V}{x_1}{\dot{x}}_1+\frac{\partial V}{x_2}{\dot{x}}_2$

$\frac{\partial V}{x_1}=a^2{x}_1+a{x}_2+x_1\Psi \left(x_1\right)$

$\frac{\partial V}{x_2}=a{x}_1+2x_2$

$\dot{V}=a\left(a{x}_1+x_2\right)+x_1\Psi \left(x_1\right)x_2+\left(a{x}_1+2x_2\right)\left(-a{x}_1-x_1\Psi \left(x_1\right)\right)$

$\dot{V}=a^2{x}_1{x}_2+a{x}_2^2+x_1{x}_2\Psi \left(x_1\right)-a^2{x}_1^2-a{x}_1^2\Psi \left(x_1\right)-2a{x}_1{x}_2-2x_1{x}_2\Psi \left(x_1\right)$

The terms listed above sum up to this;

$\dot{V}\left(x_1,x_2\right)\le -ab{x}_1^2-a{x}_2^2-b{x}_1{x}_2,$

$\dot{V}\left(x_1,x_2\right)\le -\left(ab-1\right)x_1^2-{\left(b{x}_1-x_2\right)}^2-\left(a-1\right)x_2^2$

It follows that

$\dot{V}\left(x_1,x_2\right)\le 0,$

That is $\dot{V}\left(x_1,x_2\right)=0$ if and only if $\left(x_1,x_2\right)=\left(0,0\right)$ and $\dot{V}\left(x_1,x_2\right)<0$ if $\left(x_1,x_2\right)\ne \left(0,0\right)$ . Thus, $x_1,x_2$ in the nonlinear system (2) is asymptotically stable based on Lyapunov’s theory under the imposed assumptions (i), (ii) as the time $t\to \infty$ . Since the Lyapunov function is positive definite and its derivative is negative definite, the system is asymptotically stable, meaning it will eventually return to equilibrium. Next, we introduce the bounding effect of the solutions of the nonlinear Equation (2).

5. Boundedness of Solution

The boundedness of solutions of a nonlinear differential equation of motion (2) refers to whether the solutions (e.g., the angular position, velocity, or other variables) remain within a finite range over time. In other words, it addresses whether the motion described by the differential equation will stay within certain limits and not grow indefinitely.

A solution $X\left(t,x_1,x_2\right)$ of Equation (2) is said to be bounded if exists a constant $K>0$ such that $\left|X\left(t,x_1,x_2\right)\right|<K$ , where K may depend on each solution

Theorem 1. In addition to the assumptions imposed on the parameters in Theorem 1, let

(iii) $ab>0$

(iv) $\left|\delta T\left(t,x_1,x_2\right)\right|\le \chi \left(t\right)\left(\left|x_1\right|+\left|x_2\right|\right)$

where $\chi$ is a positive continuous function of t, ${\displaystyle {\int }_0^t\chi \left(s\right)\text{d}s}\le D<\infty$ , $D>0$ , a constant, then there exists a constant $K>0$ such that every solution of (2) satisfies $\left|x_1\right|<K$ , $\left|x_2\right|<K$ , $\forall t>0$ .

Proof Also, (2) is linearized at the presence of $\delta T\left(t,x\right)$ to yield the form

${\dot{x}}_1=x_2$

${\dot{x}}_2=a{x}_2-bG\left(x_1\right)-c\sin \left(x_1\right)+\delta T\left(t,x_1,x_2\right)$

Using the Lyapunov function, V, as a scalar function for the non-homogeneous part of the system (2), that is, $\delta T\left(t,x_1,x_2\right)\ne 0$ , defined by

$2V\left(x_1,x_2\right)={\left(a{x}_1^2+x_2\right)}^2+x_2^2+4{\displaystyle {\int }_0^{x_1}s\Psi \left(s\right)\text{d}s}$ (5)

This can be re-written as

$V\left(x_1,x_2\right)\ge x_2^2+{\left(a{x}_1+x_2\right)}^2+x_1^2\Psi \left(x_1\right)$

$V\left(x_1,x_2\right)\ge x_2^2+{\left(a{x}_1+x_2\right)}^2+b{x}_1^2$

Time derivatives include

$\dot{V}=\frac{\partial V}{x_1}{\dot{x}}_1+\frac{\partial V}{x_2}{\dot{x}}_2$

$\frac{\partial V}{x_1}=a\left(a{x}_1+x_2\right)+x_1\Psi \left(x_1\right)$

$\frac{\partial V}{x_2}=a{x}_1+2x_2$

$\dot{V}=a\left(a{x}_1+x_2\right)+x_{1}G\left(x_1\right)x_2+\left(a{x}_1+2x_2\right)\left(-a{x}_1-x_1\Psi \left(x_1\right)+\delta T\left(t,x_1,x_2\right)\right)$

$\begin{array}{l}\dot{V}=a^2{x}_1{x}_2+x_2^2+x_1^2{x}_{2}G\left(x_1\right)-a^2{x}_1{x}_2-a{x}_1^2\Psi \left(x_1\right)\\ -2a{x}_2^2-2x_1{x}_2\Psi \left(x_1\right)+\left(a{x}_1+2x_2\right)\delta T\left(t,x_1,x_2\right)\end{array}$

Since

$\left|x_1{x}_2\right|\le \frac{1}{2}\left(x_1^2+x_2^2\right),\left|G\left(x_1\right)\right|\le b,\left|x_1\right|\le \left(1+x_1^2\right),$

$\left|\delta T\left(t,x_1,x_2\right)\right|\le \chi \left(t\right)\left(\left|x_1\right|+\left|x_2\right|\right)$

$\begin{array}{l}\dot{V}=a^2{x}_1{x}_2+a{x}_2^2+x_1{x}_2\Psi \left(x_1\right)-a^2{x}_1^2-a{x}_1^2\Psi \left(x_1\right)\\ -2a{x}_1{x}_2-2x_1{x}_2\Psi \left(x_1\right)+\left(a{x}_1+2x_2\right)\delta T\left(t,x_1,x_2\right)\end{array}$

$\dot{V}\left(x_1,x_2\right)\le -ab{x}_1^2-a{x}_2^2-b\left|x_1{x}_2\right|+\left(a{x}_1+2x_2\right)\delta T\left(t,x_1,x_2\right),$

$\dot{V}\left(x_1,x_2\right)\le -ab{x}_1^2-b\left|x_1{x}_2\right|-a{x}_2^2+\left(a\left|x_1\right|+2\left|x_2\right|\right)\left(\left|x_1\right|+\left|x_2\right|\right)\chi \left(t\right)$

$\dot{V}=ab{x}_1^2-\frac{b}{2}{\left(x_1^2+x_2\right)}^2-a{x}_2^2+\left(a{x}_1+2x_2\right)\delta T\left(t,x_1,x_2\right)$

$\begin{array}{l}\dot{V}\left(x_1,x_2\right)\le -ab{x}_1^2-\frac{b}{2}\left(x_1^2+x_2^2\right)-a{x}_2^2+\left(a+a{x}_1^2+2+2x_2^2\right)\chi \left(t\right)\\ +\left(2+x_1^2+x_2^2\right)\chi \left(t\right)\end{array}$

$\dot{V}\le k_2\chi \left(t\right)+k_1\left(x_1^2+x_2^2\right)\chi \left(t\right)$

where

$k_2=\max \left\{a,\left(a+2\right),2\right\}$

$\dot{V}\le k_2\chi \left(t\right)+k_1\chi \left(t\right)V\left(t\right)$ (6)

integrating the inequality from 0 to t, we have

$V\left(t\right)-V\left(0\right)\le k_2{\displaystyle {\int }_0^t\chi \left(s\right)\text{d}s}+k_1{\displaystyle {\int }_0^{t}V\left(s\right)\chi \left(s\right)\text{d}s},$

By Gronwall inequality from results of Olutimo and Omoko, ([12] [13])

$V\left(t\right)\le k_3\exp \left(k_2{\displaystyle {\int }_0^{t}V\left(s\right)\chi \left(s\right)\text{d}s}\right)=K,$

Thus, the solution of (2) satisfies

$\left|x_1\right|\le K,\left|x_2\right|\le K$

provided $\left|\Psi \left(x_1\right)\right|>b$ , $ab>0$ , K is a positive constant.

If the solution does not satisfy this condition, it is considered unbounded, meaning the variables can grow without bounds, which would typically indicate an unrealistic or unstable physical system.

6. Simulation and Results

As a work example, we are to analyze the systems of the nonlinear differential equation of motion (2). Nonetheless, we use Wolfram Mathematica and Matlab tools to solve replacing the parameters with appropriate values.

Relating our data with those present in the literature, Amélie Chevalier et al. [6], we have (see Table 1)

Table 1. Equation of motion parameter.

Parameters	Values
Shank Inertia, I	0.44
Viscous Damping, K	41.1 Nm/rad
Joint Stiffness, B	6.4 Nms/rad
Gravity, g	9.8 m/s2

From the system (2) and Amélie Chevalier et al. [5], we will extrapolate all torques to be:

• a torque brought on by gravity

• a torque brought on by the shank inertia

• a torque as a result of the effect of viscous damping

• a torque introduced by joint stiffness

• applied torque.

We will solve a type of differential equation of motion (2) of the form;

$0.04\frac{{\text{d}}^{2}x\left(t\right)}{\text{d}{t}^2}+1.75\frac{\text{d}x\left(t\right)}{\text{d}t}+41.1\frac{1}{1+x{\left(t\right)}^2}+6.4\sin \left(x\left(t\right)\right)=0.5\cos t$ (7)

7. Application of Wolfram Code

Using Wolfram Mathematica code to solve while replacing the $x_i$ with variables w, x, y z since we have to linearize the system in terms of various values of the input, damping coefficients, times, etc., we have

Although the solutions of the systems exist, we are concerned with the behavior in a variety of ways, including the vector field direction to indicate asymptotic stability, the presence and absence of viscous damping, reactions when there is an external effect, and its bounding effect.

Figure 1. Represents the phase portrait of system (2).

Figure 1 depicts a phase portrait which is spiral and is asymptotically stable. Therefore, the graphical display of this system is not a line, but a vector field where the transpose of $\left({\dot{x}}_1,{\dot{x}}_2\right)$ is plotted against $\left(x_1,x_2\right)$ .

Figure 2. Simulation of the equation of motion (2) with damping.

Figure 3. Simulation of the equation of motion (2) without damping.

Figure 4. Simulation of the equation of motion (2) with external input.

Figure 5. Boundedness of the shank knee angles at $\delta =0.1,0.25,0.5,1.2$ .

Figure 2 also illustrates that if either damping is increased or decreased, shank motion is also affected. Damping helps to stabilize the joint by opposing rapid or excessive motion. In the context of the nonlinear differential equation, damping can prevent the joint angle from deviating too far from its equilibrium position, thus maintaining stability. The effect of damping is to ensure that, after any perturbation, the joint returns to a stable state without excessive oscillations or divergence, whereas Figure 3 proves that in the absence of damping, the motion continues unabated. Nonetheless, Figure 4 explains that varying the input system presents a stable system since a periodic solution is discovered. The input effect has a substantial and steady effect on the motion of the shank. The motion of the system can also be quite small because of the small amplitude of the external force. Figure 5 shows that the angles are bounded. The bounding effect is weaker earlier at a lower amplitude and abruptly advances with a higher amplitude as time increases. This is an improved and enhanced stable condition for the shank knee motion. Stability by Lyapunov Direct Methods is vital in biomechanics to ensure that knee joint movements are safe and controlled, especially in the design of medical devices and rehabilitation strategies. Boundedness refers to the property that the solutions of the nonlinear differential equation (e.g., describing joint angles, velocities, etc.) remain within a finite range over time. In biomechanics, ensuring boundedness in joint motion is crucial for normal function, injury prevention, and rehabilitation.

8. Conclusions and Suggestions

A rudimentary model of the shank movement around the knee joint is analyzed in this work using the Lyapunov methods to determine the stability and boundedness of the solutions. By adjusting the pace at which an external force acts on the knee joint, the angle and the movement of the shank are regulated. Additionally, there is the damping effect. Damping reduces the natural frequency of the system. In the case of the knee joint, this means that the joint’s response to sudden forces or movements will be slower and more controlled. This reduction in frequency helps protect the joint from resonant oscillations, which could amplify motion and increase the risk of injury.

A controlled input system affects how well the motion surrounding the knee joint performs, potentially producing large or little amplitude. Furthermore, we demonstrated how the response produces a continuous periodic solution without dying out when the system is not dampened. In other words, controlled input systems also have a bounding effect. The Lyapunov methods are effective tools for analyzing the stability and boundedness of solutions of nonlinear models of the shank movement around the knee joint. Practically, in the real world, the Stability analysis helps in understanding whether a knee joint, after a small disturbance (like a sudden shift in body weight or a minor impact or a small applied torque), will return to a stable position or whether it might lead to instability, such as injury. Also, understanding the boundedness of the knee joint’s motion can help in designing training programs, braces, and surgical interventions that keep joint movements within safe limits. For example, a knee brace might be designed to restrict motion to a bounded range that prevents ligament strain.

Acknowledgements

Sincere and special thanks to Chinoye. G. Jones, USA, for her immense contribution. We always acknowledge the support of A.L. Olutimo (Ph.D).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	BAPO (2014) A Career Guide to Prosthetics and Orthotics. The British Association of Prosthetists and Orthotists.
[2]	Better Health Channel (2024) Sport Injuries. Department of Health, State Government of Victoria, Austrialia Singapore.
[3]	Wang, S. (2022) Biomechanical Analysis of the Human Knee Joint. Journal of Healthcare Engineering, 2022, 1-8. https://doi.org/10.1155/2022/9365362
[4]	Istabul Medical Association (2019) Health, Health Tourism Agency in Türkiye.
[5]	Amélie Chevalier, A., Ionescu, C.M. and de Keyser, R. (2014) Nonlinear Predictive Control of Shank Movement around the Knee Joint. IFAC Proceedings Volumes, 47, 9265-9270. https://doi.org/10.3182/20140824-6-za-1003.00861
[6]	Amélie Chevalier, A., Ionescu, C.M. and de Keyser, R. (2014) Nonlinear Predictive Control of Shank Movement around the Knee Joint. Proceedings of 19th World Congress 2014, Cape Town, 24-29 August 2014, 9270-9286.
[7]	Sukhman, K., Sumeet, S., Swati, S. and Bertozzi, L. (2016) Internal Model Control Based Fractional Order PID Controller for Knee Joint Motion. 2016 International Conference on Industrial and Information Systems, Roorkee, 3-4 December 2016.
[8]	Maquet, P.G.J. (1984) Biomechanics of the Knee with Application to the Pathogenesis and the Surgical Treatment of Osteoarthritis. Springer.
[9]	Winter, D.A. (2009) Biomechanics and Motor Control of Human Movement. Wiley. https://doi.org/10.1002/9780470549148
[10]	Zavatsky, A.B. (1997) A Kinematic-Freedom Analysis of a Flexed-Knee-Stance Testing Rig. Journal of Biomechanics, 30, 277-280. https://doi.org/10.1016/s0021-9290(96)00142-x
[11]	Ogundare, B.S. (2009) Qualitative and Quantitative Properties of Solutions of Ordinary Differential Equations. Ph.D Thesis, Department of Pure and Applied Mathematics University of Fort Hare Alice, South Africa.
[12]	Olutimo, L.A., Oni, O.J. and Omoko, I.D. (2020) On the Boundedness Analysis of the State Variables for a Loaded DC Servo Motor System. 2020 International Conference in Mathematics, Computer Engineering and Computer Science, Ayobo, 18-21 March 2020, 1-4. https://doi.org/10.1109/icmcecs47690.2020.240899
[13]	Olutimo, A.L. and Omoko, I.D. (2020) Stability and Boundedness Analysis of a System of RLC Circuit with Response. International Journal of Applied Mathematics, 33, 39-48.
[14]	Olutimo, A.L. and Omoko, I.D. (2020) The Problem of Convergence of Solutions of Certain Third-Order Nonlinear Delay Differential Equations. Differential Equations and Control Processes, No. 1, 12-29.
[15]	Omoko, I.D. (2023) The Qualitative Properties of Solutions of Aizermann-Type Systems of Second Order Matrix Differential Equations. PhD Thesis, the School of Post Gradute Studies, Lagos State University, Ojo.