Physiological Cybernetics: An Old-Novel Approach for Students in Biomedical Systems

1. Introduction

Wiener in a seminal book (Wiener, 1948) associated the ancient Greek word ‘κυβερνητικος’ to the control of physiological systems. “Thus, as far back as four years ago, the group of scientists about Dr. Rosenblueth and myself had already become aware of the essential unity of the set of problems centering about communication, control and statistical mechanics, whether in the machine or in the living tissue. [...] We have decided to call the entire field [...] by the name Cybernetics, which we form from the Greek κυβερνητης or steersman. In choosing this term, we wish to recognize that the first significant paper on feed-back mechanisms is an article on governors, which was published by Clerk Maxwell in 1868 and that governor is derived from a Latin corruption of κυβερνητης. We also wish to refer to the fact that the steering engines of a ship are indeed one of the earliest and best developed forms of feed-back mechanisms.”

The increasing knowledge in each sector of science led to a huge diversification of scientific research, especially in a borderline sector like cybernetics applied to physiological systems. A first problem to solve was the following: let’s suppose that two groups, one with a control engineering experience and the other one with a medical background (e.g., physicians), decide to cooperate, because they strongly believe that a joined research could be useful for developing mathematical and statistical models. Usually physicians do not have enough time to study and apply advanced modelling.

Wiener approached the communication between scientists belonging to different disciplines: “If a physiologist, who knows no mathematics, works together with a mathematician, who knows no physiology, the one will be unable to state his problem in terms that the other can manipulate, and the second will be unable to put the answers in any form that the first can understand. [...] The mathematician need not have the skill to conduct a physiological experiment, but he must have the skill to understand one, to criticize one, and to suggest one. The physiologist need not be able to prove a certain mathematical theorem, but he must be able to grasp its physiological significance and to tell the mathematician for what he should look.”

A correct interaction in terms of a clear communication and reciprocal comprehension of the objectives of the research activity between groups with different competences is a crucial aspect in any interdisciplinary research.

In 2003 at the University of Pisa it was decided to introduce a new course for undergraduate students in biomedical engineering, based on the Wiener ‘utopia’, in order to teach a novel discipline useful for helping biomedical students to communicate and cooperate effectively with physicians. We named this new course as Physiological Cybernetics, remembering the old Wiener definition.

The organization of this course was a difficult task, and it required to gain experience in order to integrate so different disciplines and to produce a common language between students in biomedical engineer and physicians. At a first glance this attempt seemed to be too ambitious, because the different approaches of biomedical engineers with respect to physicians seemed incompatible and even the languages of the two groups were so different to remember the Babel tower…

A great deal of effort and attention was required to produce appealing and stimulating lectures, but after many years we can affirm that this challenge is successful, especially for the enthusiastic answers of the students: their number was increasing year after year (about seventy students per year are now attending the course).

A strict and trusted cooperation between different groups of physicians is growing up and several groups of physicians belonging to different medical fields are going to join us for new interactions.

The aim of this chapter is to describe how the approach to physiological cybernetics has led to integrate academic lessons with research activities. To be more specific, the basic idea of Physiological Cybernetics was to search for models able to emulate physiological systems based on the feedback theory and/or the system theory.

In fact, recently, the widespread use of friendly software packages for modelling, along with the development of powerful identification and control techniques has led to a renewed interest in control (Khoo, 2011; Hoppensteadt & Peskin, 2002; Cobelli & Carson, 2008) and identification (Westwick & Kearney, 2003) of physiological systems. Unfortunately physiological systems are intrinsically time variant and highly non linear. Moreover, an effective balance of the model complexity is a difficult task: low order models are usually too simple to be useful, on the other hand high order models are too complex for simulation purposes and they have too many unknown parameters to be identified.

Each model selected for investigation was studied by a group of biomedical students supervised by physicians. Each model required several iterations and reformulations, due to the continuous adjustment of the research objectives, changing their final horizon, because of the gap between experimental data and theoretical models, so that the answers to the doubts and questions were continuously moving with the obtained partial results.

A final goal of the research was to apply a mathematical framework for helping medical diagnostic techniques and for testing new therapeutic protocols.

The procedure of model extraction followed two main pathways: the first one (pathway A) led to a formulation of a mathematical model usually based on differential equations and on an as deep as possible insight into physiological mechanisms (Marmarelis, 2004; Ottesen et al., 2004; Edelstein-Keshet, 2005; Jones et al., 2009) via a physical description of the system.

The second one (pathway B) was founded on a model description based on a black-box and data-driven identification (Westwick & Kearney, 2003; Cobelli & Carson, 2008), usually leaving to stochastic models with a parametric or non-parametric structure (Ljung, 1987), depending on the a-priori knowledge of constitutive laws governing the observed system.

In this paper we will describe two examples of research activity based on the Physiological Cybernetics, both of them addressed to produce a biomedical framework for predicting the effects of therapeutic actions, but following the two different pathways. The first example follows a statistical non parametric approach, the second one a mathematical model based on differential equations

2. Therapies in obese patients: A statistical approach

In 2004, some lessons of the Physiological Cybernetics course were addressed to describe metabolic dynamics of thyroid hormones T₃ and T₄.

It was an emblematical example regarding physiological feedback theory, intrinsically embedded inside human body.

We decided to focus some lessons on the model presented by Di Stefano et al. (1975), in an interesting paper, showing how this hormonal regulatory system could be described in terms of differential equations. This item was so intriguing for students, to require the support of physicians belonging to the Department of Endocrinology and Metabolism of the University of Pisa, in order to gain a better understanding of the physiology related to the feedback regulation of thyroid hormones. It was a representative example of pathway A, typical of classic physiological feedback, with a controller -the thyroid gland- embedded in the human body.

One of the physicians proposed a different challenging test to students: how to model another pathology with growing interest in endocrinology, i.e. the obesity?

This challenge was very complex and unsolved from a mathematical viewpoint. It was a classical example of Babel tower, because what physicians expected from us was impossible to be fulfilled in a deterministic framework, similar to the approach leading to the thyroid model.

First, we tried to consider differential equations for modelling dynamics of hormones, like leptin and ghrelin, playing an important role in controlling our weight, but the results obtained were too qualitative, simple and poor to mimic the multi-factorial aspects of obesity. It seemed to be a failed attempt, because it produced a useless model.

Hence we decided to change our approach to the challenge: if a deterministic model was inadequate, a data-driven black box model could be an alternative solution and we decided to follow pathway B. We came to the conclusion that the first and reachable step for coping with obesity was to build an interactive, user-friendly and graphically oriented toolbox for classifying obese patients. Therefore a SW tool, named Obefix, was developed for helping physicians in the classification of obese patients from physiological and psychological data.

Obefix program (Landi et al., 2007) was designed in order to produce an easy-to-use software tool for capturing all essential information on the patients using a reduced data set, solving the problem related to the high data dimensionality.

An interesting outcome was that this software tool was able to classify patients in a limited and user-selected number of clusters.

Consider to analyze a numerous group of patients. First Obefix’s user may use the toolbox for searching a blind unsupervised partition of the treated data in different clusters, using a reduced set of variables, valuable for a correct classification of the patients.

After this first step, a supervised action is possible: physicians, after an evaluation of the unsupervised classification, can ask Obefix to repeat the analysis on a restricted subset of the initial individuals, in order to eventually exclude out-of-range patients (the outliers).

In this framework, physicians can easily load data, select variables of interest, run a fast analysis and visualize results. Clusters are represented in planes, the principal planes, and single patients can be followed, automatically classified as belonging to a cluster, and grouped in Excel spreadsheets.

Obefix employs PCA (Principal Component Analysis) (Jolliffe, 2002) as an engineering statistical tool for reducing data dimensionality: users can then select either hierarchical or k-means clustering methods, for classification of patients on selected principal planes.

A clinical example of Obefix application was presented in Landi et al., (2007) the case study was the a-posteriori analysis of a dataset of severe obese women, submitted to adjustable gastric banding surgery. Obese individuals were initially candidate for gastric bariatric surgery; a presurgical preparation included also psychological evaluation.

At first, Obefix toolbox was applied for a multiple regression analysis (Mardia et al., 1979) with delta BMI (variation of the Body Mass Index expressed in %) six months after the gastric banding surgery as a dependent variable, associated with changes in pre-operative psychological data tests as independent variables.

The administrated questionnaire included 567 statements and subjects had to answer ‘‘true’’ or ‘‘false’’ according to what was predominantly true or false for them. It must be remarked that these results have been obtained using only psychological data and that in the literature the quantitative extraction of effective similarities in groups of patients in the case of a so complex and multi-factorial pathology is considered a critical and unsolved problem.

Three main homogeneous clusters were identified, representing subgroups of patients with working problems, with antisocial personality disorder and with obsessive-compulsive disorder. A strict correlation was statistically verified between the variations of BMI six months after surgery with the patients belonging to each subgroup.

All conclusions regarding the similarities between individuals belonging to different clusters were in a good accordance with medical experience and with clinical literature. Since Obefix development was considered a winning experience, we proceeded toward a following step, more interesting for the aims of the physiological cybernetics, i.e., produce and use a model able not only to classify the patients, but also to predict individual therapeutic outcome in terms of Excess Weight Loss (EWL, another common index for evaluating the loss of weight) after two years from surgery, using a set of pre-surgical data.

To be clearer, the more interesting aspect of this research was to set up a software tool able to predict the effects of a therapy and to address clinical researchers in choosing the patients that will maximally benefit from surgery.

A success in this task could represent the demonstration that the novel vision of Wiener was not a utopia, but a first example of dream coming true.

The research was again addressed to the study of the loss of weight for patients submitted to adjustable gastric banding surgery, because it was intriguing to consider a case study characterized by a high level of uncertainty in the prediction of long term effects.

Nowadays, in the medical literature it is still debated which categories of patients are better suited to this type of bariatric procedure and the selection of candidates for gastric banding surgery doesn't follows standardized guidelines.

In order to create a predictive model, the use of Artificial Neural Networks (ANNs) (Bishop, 1995; Rojas, 1996) appeared to be the best solution for predicting the weight loss after bariatric surgery, with respect to more traditional and used mathematical tools, e.g., the multiple linear regression. Therefore, a particular ANN was developed (see Figure 4) to improve the predictability of the linear model using a multi-layer Perceptron (MLP) with non linear activation functions (Rumelhart et al., 1986).

A preliminary study on the feasibility of the statistical approach for obese patients was presented in Landi et al., (2010) while, a paper considering the application of ANNs in the outcome prediction of adjustable gastric banding in obese women was published in Piaggi et al., (2010).

In the following, an outline on the engineering approach to this predictive tool is briefly sketched.

The first step was to select the most significant predictors of long term weight loss (the dependent variable) among the psychological scales, age and pre-surgical BMI (independent variables) (Van Hout et al., 2005).

In order to choose the most predictive inputs of a ANN with a limited data set and several potential predictors, a best-subset algorithm based on multiple linear regression (Neter, 1975) was employed. Namely, all combinations of the independent variables (subsets including from one to four variables, in order to avoid over-fitted solutions due to a large number of parameters, with respect to observations) were separately considered as models for computing the best linear fit of the dependent variable.

The best predictive subset was selected from all these models as that with the highest adjusted R² and a p-value less than 0.05.

The result was that age and the three psychological scales Paranoia - Pa, Antisocial practices - Asp and type-A behaviour - TpA constituted the best subset, and a predicted weight loss (WL) score was estimated through the formula

based on the linear combination of their regression coefficients, i.e., regression coefficients of (1) were a measure of the linear relationship between each independent variable and WL.

A non linear model was then built upon the same variables: the aim was to increase the goodness of prediction, taking advantage of ANNs data fitting capability.

For doing this, the four selected variables were used as inputs of a standard MLP for obtaining a non linear predictive score named u (see Figure 5).

A non linear activation function (i.e., the hyperbolic tangent function) was employed at the hidden layer units of the MLP to obtain a non linear combination of the inputs, as following:

This ANN architecture extended the regression performance of the previous linear model, which can be still obtained by replacing the nonlinear activation functions with the identity functions in the MLP, removing the nonlinear capability of the model.

The u score was then obtained as:

The global cost function - minimized by the ANN training process - was based on the correlation between u and WL scores, including their standardization terms, as following:

In this way, the ANN found the optimal values of weights (W_xh and W_hu) and bias (b_xh and b_hu), which accounted the maximum correlation between the two scores.

The non linear solution accounted for 36% of WL variance, significantly higher than 10% of the linear model using the same independent variables: this indicated a better fit for the non linear model.

Furthermore, subjects were assigned to different groups according to actual WL quartiles in order to evaluate the classification (ROC curves) and prediction (cross-validation) capabilities of the estimated models. In Figure 6, the sensitivity and specificity of both models in relation to WL outcome are plotted for each possible cut-off in the so-called ROC curves and the Area Under each ROC Curve (AUC) is estimated. AUC measures the discriminating accuracy of the model, i.e., the ability of the model to correctly classify patients in their actual quartile of WL.

As a result, the non linear model achieved better results in terms of accuracy and mis-classification rates (70% and 30% vs. 66% and 34%, respectively) than the linear model.

So far, both linear and nonlinear predictive models were built by considering all patients of the data set, i.e., each model was estimated from a database with known input and output data.

After this model-building step, the linear and nonlinear models were applied to new patients, with unknown output values, in order to have a quantitative check on the effectiveness of the proposed method on the correct selection of the therapeutic effects.

Two additional statistic tools were introduced: the cross-validation method and the confusion matrix.

Both in case of linear and nonlinear model, patients were randomly subdivided in two groups, used for building and testing the models. A training data set was considered for calculating linear regression coefficients in the case of linear model and for selecting the optimal weights and bias in the case of the MLP. A test data set was used to make a prediction of the WL two years after the bariatric surgery.

Confusion matrix was the tool used for the validation of the prediction. The cross-validation method was repeated 100 times, changing the subsets of patients for training and test sets. It was surprising to verify that after this blind test on the whole dataset, it was possible to establish with over 70% of reliability if the patients will either maximally or minimally benefit from the intervention after two years, in the case of the nonlinear model. Conversely, the reliability was reduced of about 30% in the case of the linear model (Piaggi et al., 2010). Considering that the analysis was restricted to psychological presurgical tests and to age, this result seems to be a surprising success of a research derived from the physiological cybernetics course.

3. Therapies in HIV disease: A predictive control approach

The second example shows the application of model predictive control (MPC) for an optimization of the therapy in HIV disease. It applies the subject of a group of lessons held during the physiological cybernetics course, in which the predictive control theory was presented to students as an effective tool for helping (and emulating) physicians in the selection of an optimal therapy, based on the patients' responses.

The origin of this activity was born when some students asked to study a mathematical model for HIV.

It was easy to find HIV models existing in literature: many of them are well known and accepted from mathematical and from biomedical engineers as gold standards for studies in viral models.

In the literature, (Wodarz & Nowak, 1999) the simplest model presented for mathematical modelling of HIV considers only three state variables and it is mathematically described by:

System (5) consists of three differential equations. The state variables are: x, the concentration of healthy CD4⁺ T-cells; y, the concentration of HIV-infected CD4⁺ cells; v, the concentration of free HIV copies.

Healthy cells have a production constant rate λ and a death rate d. Infected cells have a death rate a, free virions are produced by the infected cells at a rate k and u is their death rate. In the case of active HIV infection the concentration of healthy cells decreases proportionally to the product xv, with a constant rate β representing a coefficient that depends on various factors, including the velocity of penetration of virus into cells and the frequency of encounters between uninfected cells and free virus.

A five-state model was developed in Wodarz & Nowak (1999). This model offers important theoretical insights into immune control of the virus based on treatment strategies, which can be viewed as a fast subsystem of the dynamics of HIV infection. It is mathematically described by:

Two states are added to (5) to describe the dynamics of w, the concentration of precursor cytotoxic T-lymphocytes (CTLp) responsible for the development of immune memory and z, the concentration of effector cytotoxic T-lymphocytes (CTLe) responsible for killing virus-infected cells cytotoxic T-lymphocyte precursors CTLp.

In the fourth and fifth differential equations c, q, b and h are relative production rate, conversion rate to effector CTLs, death rate of precursor CTLs, and of effector CTLs, respectively.

This model can discriminate the trend of infection as a function of the rate of viral replication: if the rate is high a successful immune memory cannot establish; conversely, if the replication rate is slow, the CTL-mediated immune memory helps the patient to successfully fight the infection.

In Landi & al. (2008) model (6) was modified as:

Model (7) differs from previous W-N in the new state variable r, an index of the aggressiveness of the virus, which substitutes the constant β.

An arbitrary assumption is that r increases linearly with time in untreated HIV-infected individuals, with a growth rate that depends on the constant r₀ (a higher r₀ value indicates a higher virulence growth rate). This hypothesis was verified consistent with the simulation results obtained in the case of infected people who do not show significant disease progression for many years without treatment (long-term non Progressors -LTNP).

Different typologies of patients may require to change the law describing the aggressiveness dynamics. We evaluated the possibility to adapt the model (7) to patients with different clinical progressions, changing the values of some parameters. In particular, we supposed to vary the coefficients b and h, which represent the death rate of immune defensive cells (effector CTLs and precursor CTLs). We considered the two extreme cases for HIV progression (see Figure 7): the lower values correspond to the model dynamics of LTNP patients; the higher values model the dynamics of fast progressor patients (FP).

The coefficients μ_T and μ_P represent the drug effectiveness weights for specific external inputs f_T and f_P, which represent the drug uptakes in case of Highly Active Antiretroviral Therapy (HAART).

HAART is a combination therapy that includes:

• Reverse Transcriptase Inhibitors (RTI), to prevent cell-to-cell transmission, inhibiting reverse transcriptase activity.

• Protease Inhibitors (PI), to prevent the production of virions by infected cells, inhibiting the production of viral protein precursors.

In different models presented in literature, the effects of RTI and PI drugs have been aggregated, nevertheless we decided to mimic the effects of PI drugs reducing the rate of virus production, i.e., modifying the rate coefficient k of production of new infected cells in the dynamical equation. Instead the effect of RTI drugs is simulated by reducing the infection rate of CD4⁺ cells by free virus. So, in model (7) the RTI drugs act in virulence equation, because their main role is halting cellular infection.

Another important feature differentiating the proposed model from standard literature is that it does not admit stable steady states, since the model parameters are such that, i.e., the aggressiveness never becomes a constant, since a slow increase of r describes well a real progression of the HIV infection. This hypothesis originates from the observation that the possibility of eradicating completely the virus has not been demonstrated and the HIV disease cannot be long-term controlled.

The inclusion of aggressiveness as a new state variable represented the main outcome of the study: this simple extension to Wodarz & Nowak models allowed us to mirror the natural history of HIV infection and to introduce a new state equation useful for introducing in the model the effects of pharmacologic control.

In Fig. 8 are shown the time courses of CD4 cells and virions obtained in simulation with model (7); for a qualitative validation of the model, compare the results with the plotted experimental data shown in Fig. 9 (Abbas et al., 2000).

A straightforward application of the control theory to model (7) was proposed in Pannocchia et al., (2010), with the application of a MPC strategy in anti-HIV therapy.

MPC algorithms (Mayne et al., 2000) utilize a mathematical model of the system to be controlled, to generate the predicted values of the future response. Predicted values are then used to compute a control sequence over a finite prediction horizon, in order to optimize the future behaviour of the controlled system. The control sequence is chosen minimizing a suitable cost function, including a measure of the deviation of the future state variables from reference target values and a measure of the control effort, while respecting state and control constraints. In plain words, the core of the control algorithm is an optimization algorithm, keeping the controlled variables close to their targets and within suitable constraints. The first output in the optimal sequence of control actions is then injected into the system, and the computation is repeated at subsequent control intervals.

The problem was how to adapt MPC to determine the optimal drug scheduling in anti-HIV therapy.

Some examples of MPC applied to biomedical applications like control of the glucose–insulin system in diabetics (Parker et al., 1999), anaesthesia (Ionescu et al., 2008), and HIV (Zurakowski & Teel., 2006) have been presented in literature, but all applications were considered for models admitting a steady-state stable equilibrium. On the other hand, MPC emerged as the more suitable solution for solving the drug optimal administration problem in anti-HIV therapy, even if the model was unstable. MPC algorithm pursued the following logic:

1. future outputs of the control algorithm are generated by the HIV model; measurements on individual patient are considered and compared with the predictions of the model.

2. the cost function to be minimized keeps the controlled variables e.g., CD4⁺ cells and free virions concentration close to the targets and respecting suitable soft constraints on the manipulated variables.

3. the cost function of the future control movements is minimized using a sequence of future PI and RTI drugs over the chosen control horizon, but only the first element of the suggested control sequence is applied to the system.

4. at the successive decision time, the algorithm is solved again if measurements of CD4⁺ cells and free virions concentration are available and the drug sequence is updated, repeating step c)

Some practical issues were considered (see Pannocchia et al., (2010) for a detailed study), because MPC was applied considering the two different cases of continuous applications of drugs, or of a structured interruption of therapy (STI) for patients. STI is a treatment strategy in HIV-infected patients, involves interrupting HAART in controlled clinical settings, for a specified duration of time. The possible explanation of the effectiveness of this clinical protocol was an induced autovaccination in the patients. The use of STI is currently debated between clinical researchers and most studies agree that STI may be successful if therapy is initiated early in HIV infection, but unsuccessful for people who started therapy later.

Furthermore, a discrete dosage approach required to modify the control algorithm using a non linear MPC: this was due to the clinical request to maintain a maximum dosage of drugs, as in standard HAART protocol, in order to reduce the risks of virus mutations.

Some comments are mandatory to stress the results of this model based on a differential equation deterministic approach. From the viewpoint of a model builder, two different situations have to be usually considered: basal and pathological conditions. In the case of infections, like HIV, the mathematical model have to mirror the natural evolution of HIV infection, and the pathological model must be more accurate, because today it is the only one that can be validated by experimental data, since patients are all maintained under therapy. The impact of therapy into HIV models must be introduced in a way as simple as possible, if we have to satisfy the task to formulate a model suitable for use in feedback control.

Simulation results were coherent with the medical findings: the comments of clinical researchers expert in HIV therapies were essential in testing the model and for evaluating the effectiveness of the proposed control methods.

Obtaining reliable models is relevant from a diagnostic and prognostic point of view, because it allows the physician to prove the therapeutic action using the model for testing the treatment in terms of optimal dosage and administration of drugs.

In 2008, the FDA approved an in silico model of diabetes as a pre-clinical testing tool for closed loop research at the seven JDRF Artificial Pancreas Consortium sites. The overall goal of the Artificial Pancreas project was to accelerate the development, regulatory approval, health insurance coverage, and clinical acceptance of continuous glucose monitoring and artificial pancreas technology (Juvenile Diabetes Research Foundation, 2008).

We strongly believe that also a simple but reliable in silico model of HIV can lead to an acceleration of the experimental tests for a clinical acceptance of new drugs in HIV disease.

Future activity will be devoted to develop models of HIV infection, able to include the issues of drug resistance and viral mutation, key issues for the HIV studies, and the interest of many clinical researchers in our work is encouraging in the upcoming research.

4. Conclusion

The Physiological Cybernetics course represents an example of integration between different disciplines, in order to produce a common language between students in biomedical engineer and physicians. It offers students an opportunity to verify in practice how to move theoretical lectures, based on the development of mathematical models, to a practical interaction with physicians. This fact seems obvious from an educational viewpoint, but it isn't so usual in practice, because it requires a preliminary long period for preparing a common language between researchers in different fields. Judging from the students’ excellent results, if compared to students attending under-graduated courses in previous years, the example proposed was very successful.

In this chapter we presented two examples of research applications, derived from this educational experience, demonstrating that the old-novel vision of Wiener was not a utopia, and that a synergic cooperation between biomedical engineers and physicians can lead to interesting results.