Physiological reflexes and control systems

Physiological reflexes are the basis of the organism’s survival. These reflexes allow, in a completely involuntary way, to maintain homeostasis.

Variations in environmental variables or the appearance of a new signal induce a reflex response from the integration centers, aimed at restoring balance or avoiding dangerous situations.

In the following report, the typical techniques of automatic controls are applied to the analysis of the muscle strain reflex and the pupillary reflex, analyzing their behavior and stability. The stretch reflex is induced by excessive stretching of the muscle and thus an autonomous contraction is induced. Pupillary reflex predicts a change in pupil diameter in response to a change in ambient light.

These reflexes are characterized by different systems and pathways but both can go into instability if, for different pathologies, the gain of the system increases excessively.

The following report is an extract from a project carried out for the course of Modeling and Simulation of Physiological Systems, held for Medical Engineering at the University of Rome Tor Vergata.

Introduction

A physiological reflex is an involuntary response to a stimulus, mediated by elements of the nervous system. Generally, the purpose of physiological reflexes is to maintain the homeostasis of the organism and avoid dangerous situations, keeping the organism safe.

Generally, to speak of reflexes, we start from the stimulus characterizing it to the change of an environmental variable or to the variation of a particular signal. This stimulus is felt by a sensor, in physiology a sensory receptor, which monitors the internal and external environment of the organism.

The afferent pathway starts from the receptor towards the integration centers where the stimulus is processed, in a more or less complex way, processing a response that, through an efferent path, will be transmitted to the effector.

Among the physiological reflexes of the human body, the muscle strain reflex and the pupillary reflex stand out. These reflexes differ in the inducing stimulus but are united by having as a response the contraction or relaxation of muscle fibers.

These reflexes can be studied by modeling and analyzing them with automatic control theories to observe how variations in the gain of the reflex, for example due to highly debilitating pathologies, lead to instability in both systems.

The following report proposes two different models and their implementation in the Simulink environment, analyzing their behavior.

Background

A control system is an apparatus that allows you to manipulate the behavior of the control system in relation to a temporal evolution of the input quantities. That is, a set of interconnected components that work to achieve a certain response even in the presence of external stimuli [1].

We talk about regulation, or simple control, when the output is required to remain constant as the input changes. Where possible, closed loop control systems are used where the feedback block is also identified, which generates a signal proportional to the value of the controlled quantity, returning it to the controller input. Within biological organisms there are various control systems, more or less complex, some local and others widespread throughout the organism, aimed at maintaining the homeostasis of the organism.

The term homeostasis was introduced in 1865 by Claude Bernard to represent how the organism tried to keep the internal environment of the body constant [2]. The term was then coined by Walter Bradford Cannon, expanding its meaning and introducing the fact that homeostasis is not a stationary static state but is a continuous evolution of states of equilibrium and conditions that can continuously vary [3]. This is due to the high complexity of the physiological control systems that certainly involve the nervous system but also the individual organs.

Therefore, the extension of the concept of homeostasis can be seen as the oscillation over time around an average value characteristic of each physiological variable. This oscillation over time can refer to seconds, minutes, days, hours or years to seconds of which parameter is considered. However, we spend most of our time oscillating between a maximum and minimum value, around the average value, in a range considered normal. This range is precisely what is called homeostasis [4].

Neuromuscular reflex

When a whole muscle is passively stretched the neuromuscular spindle fibers stretch and increase the firing rate in the afferent nerve fibers. Their sensory endings, which terminate on the stretched intrafusal fibers, are such that the afferent neuron establishes synaptic contacts directly on the motor neuron 𝛼 which innervates the extrafusal fibers of the same muscle. The activation of this synaptic contact causes the contraction of the muscle that was excessively stretched.

This myotatic reflex acts as a local negative feedback mechanism that allows to oppose any passive variation of the muscle length, so as to maintain the optimal stretch length. Examples are the patellar reflex or the stretch reflex in the arm, both reflexes are also used as a routine clinical test to evaluate the nervous system. These reflexes may be absent or depressed due to loss of excitatory input from higher levels or due to loss of inhibitory input to motor neurons from higher brain levels [5].

This reflexes involves the muscle itself, the neuromuscular spindle and motor neurons. Such participants can be reviewed as a control system. By focusing on the upper limb, it is possible to see the reflex center in the spinal cord, where the synapse between the afferent pathways and the motor neurons takes place 𝛼 to regulate the muscle according to the disturbance induced by the load.

In order to be able to analyze the individual systems, it is necessary to adopt a descriptive mathematical model. Starting from the geometric description in fig. 2 the system is described by the following equation:

\begin{equation}
M_{x}(t)-M(t)=J \ddot{\theta}
\end{equation}

Muscle model

Hill’s muscle model is considered [6]. This model, with lumped parameters, considers the muscle as a series of two zones. The first zone is given by a parallel of a driving torque (muscular force) and a damping, and is in turn in series with an elastic stiffness (fig. 1a). While the angular displacement will be the sum of the two blocks, the resulting torque will be equal in parallel:

\begin{equation}
\begin{gathered}
M(t)=K\left(\theta-\theta_{1}\right) \\
M(t)=M_{0}(t)+B \dot{\theta}_{1}
\end{gathered}
\end{equation}

Thus, deriving the first and substituting 𝜃1 in the second, it can be dynamically described by the following equation:

\begin{equation}
M_{x}(t)-M_{0}(t)=\frac{B J}{K} \dot {\ddot{\theta}}(t)+J \ddot{\theta}(t)+B \dot{\theta}(t)
\end{equation}

And, in the domain of Laplace, it is:

\begin{equation}
M_{x}(s)-M_{0}(s)=s\left(\frac{B J}{k} s^{2}+J s+B\right) \theta(s)
\end{equation}

Neuromuscular spindle

For the muscle spindle, the Soechting model is adopted [7] where the spindle is differentiated between the internal part of the nucleus and the polar part. Consider the series of a zone given by the parallel of an elastic stiffness and a viscous damper with a second elastic stiffness (Fig. 1c). The variable of interest is always the angular coordinate and the following relations hold:

\begin{equation}
\begin{aligned}
&M_{s}(t)=K_{s s}\left(\theta-\theta_{2}\right) \\
&M_{s}(t)=K_{s p} \theta_{2}+B_{s} \dot{\theta}_{2}
\end{aligned}
\end{equation}

Reflexes

As a result, we know that what happens at the spindle level counts, that is, we will have the response linked to elongation through a certain gain:

\begin{equation}
M_{0}(t)=\beta\left(\theta-\theta_{2}\right)
\end{equation}

From which the system equation is found:

\begin{equation}
\dot{M}_{0}(t)+\frac{M_{0}(t)}{\tau}=\beta\left(\dot{\theta}\left(t-T_{D}\right)+\frac{\theta\left(t-T_{D}\right)}{\tau \eta}\right)
\end{equation}

Where the two parameters, functions of the lumped constants, are:

\begin{equation}
\tau=\frac{B_{s}}{K_{s s}+K_{s p}} \quad \text { e } \eta=\frac{K_{s s}+K_{s p}}{K_{s p}}
\end{equation}

Passing in the Laplace domain we can derive the couple acting as:

\begin{equation}
M_{o}(s)=\beta\left(\frac{s \tau+s / \eta}{s \tau+s}\right) e^{-s T_{D}} \theta(s)
\end{equation}

For the implementation in Simulink, physiological parameters known in the literature are considered [8, 9], in Tab. 1. This system can be represented in a block diagram, typical of automatic controls, shown in fig. 3.

Parameter	Value
B	2 [Nms]
J	0.1 [kg m^2]
Td	0.02 [s]
beta	100
eta	5
k	50 [Nm]
tau	1 / 300 [s]

Pupillary light reflexes

Another reflex of great clinical interest is the pupillary reflex, also known as the photomotor reflex. This reflex controls the diameter of the pupil in response to the intensity of the light reaching the retina. The optic nerve constitutes the afferent pathway as it perceives incoming light through the ganglion cells of the retina. The information reaches the pretectal nucleus in the superior midbrain and then the Edinger-Westphal nucleus. From here the oculomotor nerves branch off to the ciliary ganglia and the inversion of the constrictor muscle of the pupil [10].

Under physiological conditions, the pupils both respond identically. The comparison of the asymmetric responses allows to evaluate a lesion. For example, a direct response in one pupil without a consensual response in the other indicates a possible problem in the motor connection of the second. Or, it is possible to have direct reflexes losses while the consensual response remains active.

Stark model

For the implementation the Stark model [11] is considered.

This model interprets the flow of light as the light intensity acting on the pupil area and the controller’s goal (pupillary reflex) is precisely to maintain this constant flow. If the flow varies, the nervous system is activated to vary the area of the pupil. This model contains many non-linearities but, by linerizing it around the working point, it allows to adequately describe the physiological response [8]. The block diagram of the system is shown in fig. 4. We therefore consider an LTI system with a third-order transfer function in series with a delay with time constant 𝜏. This system was described by Stark following an accurate fitting of the experimental data.

Or the loop transfer function:

\begin{equation}
H(s)=\frac{K e^{-s D}}{(1+\tau s)^{3}}
\end{equation}

Where the gain is given by the product of the gain of the transfer function of the area change and the linearized factor of light intensity (with respect to which the model is linearized), namely:

\begin{equation}
K=I_{\mathrm{ref}} K_{1}
\end{equation}

The system parameters are present in Tab. 2.

Stability

For a dynamic system, an equilibrium state is a state in which the system remains indefinitely, in the absence of perturbations. If this state of equilibrium is stable (asymptotically), this allows the system to actually remain there.

For a linear system we can speak of stability of the system itself (and not only of the equilibrium point) and in particular it is stable if the eigenvalues of the dynamic matrix of the system are a negative real part [12]. This mathematical definition can be particularized in different criteria to analyze the stability of the single system. Considering a generic characteristic equation feedback system:

\begin{equation}
1+G(s) H(s)=0
\end{equation}

Where 𝐺 (𝑠) is the direct transfer function and 𝐻 (𝑠) the feedback transfer function and 𝐺 (𝑠) 𝐻 (𝑠) the loop transfer function (fig. 6). So such that:

\begin{equation}
Y(1+G H)=G X
\end{equation}

And then the transfer function:

\begin{equation}
W=\frac{Y}{X}=\frac{G}{1+G H}
\end{equation}

Considering that 𝐺 (𝑠) 𝐻 (𝑠) is a rational function divided by gain 𝐾 (positive for negative feedback systems), the roots of this equation describe a curve in the plane 𝑠 which is called the locus of the roots.

\begin{equation}
G(s) H(s)=K \frac{N(s)}{D(s)}
\end{equation}

To have stability, all the poles of the transfer function are required to be on the left side of the plane, ie apart from the real negative part [1]. These poles correspond to the zeros of eq. (12) or to study the trend:

\begin{equation}
D(s)+K N(s)=0
\end{equation}

Simulink

Simulink [13] is software developed by MathWorks that provides a graphical approach based on a modeling environment that allows the user to convert the problem into a network of mathematical function blocks. In addition, this environment allows integration with the Matlab language and its programming functions.

Through Simulink it is possible to implement the considered models by adding various features, such as dashboards to interact with the system parameters, useful for analyzing their behavior and stability.

Results

These models were implemented within Simulink by analyzing the stability and the influence of the different coefficients, in particular of the gain.

Neuromuscular reflexes

The Simulink model of the stretch reflex is shown in fig. 7 where, in addition to the model already described in the block diagram in fig. 3 there are some graphic blocks to change the variables involved and to save the data. In particular, the fnc block is described in the appendix (section 4.1).

The system is analyzed over a time from 0 to 1 s by observing the step response after 0.1 s. The trend of the signal is shown in fig. 8. The system moves towards the state of equilibrium remaining stable and constant. Initially there is a slight oscillation around the equilibrium value, due precisely to the gain of the system which tends to compensate for the suddenly added load.

This oscillatory effect increases as the gain increases, as indicated in fig. 10a. In particular, for a lower gain (𝛽 = 50) the system stops oscillating but tends to balance itself at a much greater angle, up to about 30 ° with respect to the initial configuration which we remember being with the elbow not strictly perpendicular but with an angle of about 135 °.

Underdamping

This under-damping is not present for higher gains where the angle of equilibrium becomes much lower but the response begins to become oscillatory. In particular, for a gain of 250 the response begins to diverge and the system is unstable.

From the stability analysis, the limit gain is identified at 202.62. This value is identified by considering the Padè expansion for the exponential and analyzing the place of the roots (fig. 11a) or using the Matlab margin() command which returns the Bode diagram and the limit gain (fig. 11b).

Finally, a further analysis (fig. 10b) shows how as the response delay increases, the value of the limit gain decreases. Therefore, by setting 𝛽 between the different simulations at the default value (𝛽 = 100), the system will be unstable for delays greater than 0.1.

Pupillary reflexes

To analyze the response of the pupillary reflex, a sine wave of frequency 𝑓 = 2𝜋 ∗ 0.7 [rad] = 0.7 Hz [11] is considered. The analysis is conducted from 0 to 10 s, analyzing the response in terms of pupillary dilation.

The Simulik model is present in fig. 12. In addition to the previous block diagram illustrated (fig. 5) there are two inputs, sinusoidal and step, used individually (commenting on one or the other), the fcn block analogous to the case of the stretch reflex and described in section 4.1, a scope to display trends and a block to save the data.

A sinusoidal input corresponds to a similar response which, for sufficiently low gain values, follows the input with a compliant, albeit out of phase, oscillation. This oscillation is affected by a distortion effect for gain values close to the limit value until it becomes unstable and diverging for values above the limit value. In particular, in fig. 15a, 5 different values of 𝐾 and the responses resulting from the same input are compared.

Frequency response

Furthermore, it is possible to observe how the unstable oscillation is at a much higher frequency than that of the input. The gain limit value is equal to 𝐾𝑙𝑖𝑚 = 1.852 and was identified with a stability analysis through the place of the roots in fig. 14.

In addition, a series of step input analyzes were also carried out to better analyze the system response and its oscillatory nature. Neglecting the default value of the gain, for which the response is stable and free of oscillations, for the values analyzed in fig. 15b, the answer fluctuates. This oscillation is damped or divergent depending on whether the gain is lower or higher than the limit value.

Oscillations

Analyzing the oscillation frequency (∈ [0.8; 1.1] Hz) this is similar to the values typically indicated in the literature for the hippus phenomenon [14] and also to the critical frequency of the system (7.13 rad / s). Hippus is a particular phenomenon of rhythmic movement of pupillary dilation and contraction that occurs in pathological conditions such as acute meningitis, epilepsy or other brain diseases. It can therefore be seen as an excessive gain of the system, such as to bring it into instability.

Conclusions

To subject the physiological mechanisms and reflexes to an analysis according to the classical criteria of the theory of automatic controls, it is necessary to obtain a model of a dynamic system that is sufficiently accurate and at the same time sufficiently simple in terms of computational costs and parameters to be estimated with respect to the organism of interest.

The implementation of a physiological reflexes controller model within an interactive simulation environment, such as Simulink, allows you to analyze it and understand the effect of the individual parameters.

Furthermore, it is possible to exploit the classical criteria for stability analysis and verify which parameters lead to instability in the system by simulating different pathological conditions.

Code availability

The code is available at the project’s GitHub repository.

References

• [1] Giovanni Marro. Controlli automatici. Italian. Bologna: Zanichelli, 2004
• [2] Claude Bernard. An introduction to the study of experimental medicine. Vol. 400. Courier Corporation, 1957.
• [3] Walter Bradford Cannon. “The wisdom of the body”. In: (1939).
• [4] Kelvin J.A. Davies. “Adaptive homeostasis”. en. In: Molecular Aspects of Medicine 49
• [5] Lauralee Sherwood. Fisiologia umana: dalle cellule ai sistemi. Italian. OCLC: 849240778. Bologna: Zanichelli, 2008.
• [6] HIll. “The heat of shortening and the dynamic constants of muscle”. en. In: Proceedings of the Royal Society of London. Series B – Biological Sciences 126.843 (Oct. 1938), pp. 136–195.
• [7] R. E. Mains and J. F. Soechting. “A Model for the Neuromuscular Response to Sudden Disturbances”. en. In: Journal of Dynamic Systems, Measurement, and Control 93.4 (Dec. 1971), pp. 247–251.
• [8] Michael C. K. Khoo. Physiological control systems: analysis, simulation, and estimation. en. Second editon. IEEE Press series in biomedical engineering. Piscataway, NJ : Hoboken, New Jersey: IEEE Press ; Wiley, 2018.
• [9] J. F. Soechting et al. “Evaluation of Neuromuscular Parameters Describing Human Reflex Motion”. en. In: Journal of Dynamic Systems, Measurement, and Control 93.4 (Dec. 1971), pp. 221
• [10] Eric R. Kandel et al., eds. Principles of neural science. Sixth edition. New York: McGraw Hill, 2021.[11] Lawrence Stark and Philip M. Sherman. “A SERVOANALYTIC STUDY OF CONSENSUAL PUPIL REFLEX TO LIGHT”. en. In: Journal of Neurophysiology 20.1 (Jan. 1957), pp. 17–26.
• [12] Osvaldo Maria Grasselli, Laura Menini, and Sergio Galeani. Sistemi dinamici: introduzione all’analisi e primi strumenti di controllo. Italian. OCLC: 898636256. Milano: Hoepli, 2008. ISBN: 978-88-2033703-2.
• [13] MathWorks. “Simulink”. In: (2022). URL: https:// it.mathworks.com/products/simulink.html.
• [14] Philip R. K. Turnbull et al. “Origins of Pupillary Hippus in the Autonomic Nervous System”. en. In: Investigative Opthalmology & Visual Science 58.1 (Jan. 2017)

Appendix

Below is a description of some important parts for the implementation of the model and the production of the graphs.

fcn block

Both reflexes models using the fcn block to automate some features. In particular, this function block calculates the limit gain and compares it with the current gain (set via the dashboards).

Depending on whether it is greater or less, it colors the display red or green. To be able to do this, unlike classic Matlab, it is necessary to declare the functions as extrinsic and initialize the new variables, also declaring their class. In addition, the block also allows you to adjust the ranges of the colored slider by dividing it into three regions (green, orange and red) and update it with each simulation.

function [calcGain,u]=fcn(u,B,J,Td,eta,k,tau)
% this avoid problem with function and code generation
coder.extrinsic(’pade’); coder.extrinsic(’tf’); coder.extrinsic(’rlocus’)
coder.extrinsic(’tf’); coder.extrinsic(’margin’); coder.extrinsic(’type’) coder.extrinsic(’strcat’); coder.extrinsic(’get_param’); coder.extrinsic(’set_param’)
% setup color as strings
green=’[ 0.2471 0.7804 0.2392]’; red=’[0.7255 0.2941 0.2100]’;
% set trasfer function
G_num =1;
G_den=[B*J/k J B]; int =[1 0];
F_num =[ tau 1/ eta ]; F_den =[ tau 1];
% calc pade for numerator exponential expression
[expon_num_pade,expon_den_pade] = pade(Td,4);
% calc transfer function
num=conv(expon_num_pade ,conv(G_num ,F_num));
den=conv(expon_den_pade,conv(int,conv(G_den, F_den))); sys=tf(num,den);
% calc gain margin
% new varaibles need to be initialized also with class type
Gm_pade=double(0); Gm_pade=margin(sys); if u>Gm_pade
set_param(’stretch_reflex/Display1’,’BackgroundColor’,red); set_param(’stretch_reflex/Display’,’BackgroundColor’,red);
else
set_param(’stretch_reflex/Display’,’BackgroundColor’,green);
set_param(’stretch_reflex/Display1’,’BackgroundColor’,green); end
% set parameters for the Linear Gauge % need struct array like 3x1
% range.Min with 3x1 minimum, range.Max with 3x1 maximum, range.Color with [r g b] color
range=struct(’Min’,{0; 0.9*Gm_pade; 1.1*Gm_pade},’Max’,{0.9*Gm_pade; 1.1*Gm_pade; 250},’Color’
,{[0.2471 0.7804 0.2392];[0.9294 0.6941 0.1255]; [0.7255 0.2941 0.2000]}); set_param(’stretch_reflex/Linear Gauge’,’ScaleColors’,range)
calcGain=Gm_pade;
end

Output

That block also acts as pass through for the variable u or the current gain so that this can be saved, at the end of the simulation, in the results block, without having to wait for the variables to be updated (which would happen at the next simulation).

The output block is renamed from simulation to simulation through a callback function, specific to the model, which is activated every time the simulation is initialized. By means of a callback function, the file is saved in the folder ../results/ and that the name contains the concatenation of the parameter name and its value, for example:

\text { 'results } /+\text { gain }+\{\text { gain }\}+\mathrm{Td}+\{\mathrm{Td}\}+\text { mat' }^{\prime}

The results are saved as Timeseries.

Post processing

Within the results there are, in addition to the input and the response, also the values of the parameters. These are used for the legends and titles of the graphs by uploading, file by file, within a given structure{} with the following fields:

file_struct=dir(’results/*.mat’);
% load all the files in data{} structure with different fields
VAR ={};
for i=1:length(file_struct) 
load(strcat(’results/’,file_struct(i).name))
data{i}.Time=stretch_reflex_result.Time;
% use mode() function to avoid numerical error
% this data have to be consant)
data{i}.Gain=mode(stretch_reflex_result.Data(:,1));
data{i}.input=stretch_reflex_result.Data(:,2); data{i}.response=stretch_reflex_result.Data(:,3);
data{i}.Td=mode(stretch_reflex_result.Data(:,4)); data{i}.margin=mode(stretch_reflex_result.Data(:,5));
% save the name of the file to export corresponding figs
data{i}.Name=erase(file_struct(i).name,’.mat’);
clear stretch_reflex_result 
end