Fractional-Order Model Predictive Control for Neurophysiological Cyber-Physical Systems: A Case Study using Transcranial Magnetic Stimulation

Fractional-Order Model Pr edictiv e Contr ol f or Neur oph ysiological Cyber -Physical Systems: A Case Study using T ranscranial Magnetic Stimulation Orlando Romero †∗ Sarthak Chatterjee ‡∗ S ´ ergio Pequito † Abstract — Fractional-order dynamical systems are used to describe processes that exhibit temporal long-term memory and power -law dependence of trajectories. There has been evidence that complex neurophysiological signals like elec- troencephalogram (EEG) can be modeled by fractional-order systems. In this work, we propose a model-based approach for closed-loop T ranscranial Magnetic Stimulation (TMS) to regulate brain activity through EEG data. Mor e precisely , we propose a model predictive control (MPC) approach with an underlying fractional-order system (FOS) predicti ve model. Furthermore, MPC offers, by design, an additional layer of rob ustness to compensate f or system-model mismatch, which the mor e traditional strategies lack. T o establish the potential of our framework, we focus on epileptic seizure mitigation by computational simulation of our proposed strategy upon seizure-lik e events. W e conclude by empirically analyzing the effectiveness of our method, and compare it with event-trigger ed open-loop strategies. I . I N T RO D U C T I O N In the context of neurophysiological signals, tempo- ral fractional properties in both health and disease states hav e become apparent and, with it, fractional-order sys- tems hav e demonstrated a huge potential for clinical applications [1]–[7]. Practically , this leads signals to become non-stationary and to possess long-term memory dependen- cies with themselves, with the backwards-decaying weights of such dependencies follo wing a po wer-la w distribution [8]. The persisting temporal dependencies illustrated by these systems have giv en rise to fr actional-or der based modeling, design, and analysis of novel neurotechnologies. Recently , dynamical spatiotemporal fractional models hav e been proposed as a tool to model neurophysiolog- ical signals suitable to deal with structured data and to equip us with modeling capabilities that capture spatial (i.e., the contributions of the signal’ s components into each other) and temporal long-range memory through the so-called fractional-order coef ficients associated with the power -law exponents [9]–[15], and possibly under unkno wn unknowns [16], [17]. Notwithstanding the abov e, the main advent of model-based approaches is that we can understand † Orlando Romero and S ´ ergio Pequito are with the Department of Indus- trial and Systems Engineering, Rensselaer Polytechnic Institute, Troy , NY 12180, USA { rodrio2,goncas } @rpi.edu ‡ Sarthak Chatterjee is with the department of Electrical, Computer , and Systems Engineering, Rensselaer Polytechnic Institute, Tro y , NY 12180, USA chatts3@rpi.edu ∗ Equal contribution how an external signal or stimulus would craft the dynamics of the process. Simply speaking, it enables us to design a sequence of interactions (i.e., a control strategy or law) with the system such that we can steer its dynamics to wards satisfying desirable properties. That said, due to the highly dynamic nature of neurophysiological processes, it is imperativ e that we consider feedbac k mechanisms [18]. In other words, we need to leverage the continuous flow of measurements of the system to tune (for the individual’ s process) the control strategy . A particularly successful strategy that has achiev ed remarkable success in sev eral engineering applications is the strategy of model pr edictive contr ol (MPC) that consists of three key ideas [19]–[24]: ( i ) a model-based approach; ( ii ) capability of predicting the ev olution of the system and its states upon a devised feedback control strategy that aims to optimize an objecti ve that encapsulates the risk assessment of abnormal behavior; and ( iii ) receding finite-horizon re-ev aluation of the control strategy performance devised in the previous point. In the conte xt of neurophysiological processes, we propose to le verage fractional-order models to equip us with the aforementioned prediction and control capabilities that go hand in hand with the closed-loop design of Cyber-Physical Systems (CPSs). As a consequence, we will be able to dev elop stimulation strategies in the form of T ranscranial Magnetic Stimulation (TMS) from electroencephalographic (EEG) data, and use it to annihilate or mitigate the o verall du- ration and strength of an epileptic seizure. Notwithstanding, we belie ve that similar design and strategies can be en vi- sioned in other conte xts where closed-loop deep-brain elec- trical stimulation is a v ailable, e.g., Parkinson’ s disease [25], [26], Alzheimer’ s disease [27], depression [28], [29], and anxiety [30], just to mention a fe w . T o summarize, in what follows we introduce in a ped- agogical manner the control mechanisms to be deployed as part of future neurophysiological c yber-physical systems, with particular emphasis on TMS for epilepsy . Ultimately , the integration of these design features will lead to more reliable CPSs that will immediately improve, or otherwise bring a positive impact, on the quality of life of the patients that qualify for the use of such technologies. The paper is organized as follows. Section II introduces the preliminary setup of our problem. Section III presents the fractional-order model predictiv e control framework. Finally , we present illustrativ e examples demonstrating the ef ficacy of epileptic seizure mitigation using TMS in Section IV. I I . P R E L I M I NA R I E S There are primarily two kinds of brain stimulation strate- gies that exist, which include open-loop and closed-loop strategies. Open-loop stimulation strategies consist of any stimulation strategy that does not utilize current brain acti vity data to regulate the stimuli applied to the patient’ s brain. On the other hand, closed-loop (brain-responsi ve) strategies con- sist of stimulation treatments based on automatic electrical stimulation directly influenced by the present (i.e., real-time) behavior being observed through continuous recording of brain activity . Past data can also be used in this strategy , but due to limited storage capabilities, these mechanisms are usually designed to depend e xclusi vely on the actual recorded and stored data at a gi ven time, consisting of a finite temporal window ranging from a fixed number of past instances of discretized time to the present measurement. A. T ranscranial Magnetic Stimulation T ranscranial Magnetic Stimulation (TMS) is a noninv asiv e form of brain stimulation created by inducing electric cur- rents at specific areas of the brain by the help of a changing magnetic field using electromagnetic induction. This is done using an electric pulse generator connected to a magnetic coil, which in turn, is connected to the scalp. Although research on this form of neurostimulation is still evolving, it has been sho wn to demonstrate therapeutic potential in neurodegenerati ve disorders like Alzheimer’ s disease [31], motor neuron disease [32], and stroke [33]. B. Linear fractional-or der systems For many biological systems, linear time-in variant (L TI) state-space models are insufficient to accurately capture the real e volution of the systems for anything other than a v ery small interval of time into the future, giv en that the current state of the system may have a non-negligible dependence on sev eral past states, or ev en from the states ranging from the entire period of time so far . In this paper , we focus solely on discrete-time systems and control, and we will model the abov e scenario as x k +1 = f k ( x k , x k − 1 , . . . , x 0 ) + w k , (1) with known functions f k : R n ( k +1) → R n for k ∈ Z + and w k being the pr ocess noise . Alternativ ely , we may consider finite-history models of the form x k +1 = f k ( x k , . . . , x k − p +1 ) + w k , (2) such as multi variate autore gressi ve (MV AR) models x k +1 = p − 1 X j =0 A j x k − j + w k , (3) with A 0 , . . . , A p − 1 ∈ R n × n . Ho we ver , system identification for such autoregressiv e models can become mathematically intractable or suffer issues of numerical instability that originate from the possibly large number of parameters to be estimated ( n 2 p entries in A 0 , . . . , A p − 1 , plus usually a certain other number for the covariance matrix of the noise). For the reasons covered above, we introduce the so-called linear fractional-or der system (FOS) models, of the form ∆ α x k +1 = Ax k + w k , (4) where the state coupling matrix A ∈ R n × n and the vector of fractional-or der exponents α ∈ R n + are now the only parameters to be estimated. This model can be rewritten [34] as x k +1 = k X j =0 A j x k − j + w k , (5) which can be readily approximated as an MV AR model akin to (3). Furthermore, for the simulations conducted in this paper , we will model w k as an additive white Gaussian noise (A WGN). I I I . F R AC T I O N A L - O R D E R M O D E L P R E D I C T I V E C O N T R O L W e start by showing that (3) can be written as an aug- mented L TI system model. T o do this, let ˜ x k =    x k . . . x k − p +1    (6) denote the so-called augmented state vector , with the un- derstanding that x k = 0 for k ∈ {− 1 , . . . , − p + 1 } . Then, clearly , the first block in ˜ x k +1 can be expressed as a linear combination of the block in ˜ x k through (3). On the other hand, the remaining p − 1 blocks in ˜ x k +1 precisely match the first p − 1 blocks of ˜ x k . More precisely , we ha ve ˜ x k +1 =      A 0 . . . A p − 2 A p − 1 I . . . 0 0 . . . . . . . . . . . . 0 . . . I 0      | {z } = ˜ A ˜ x k +      I 0 . . . 0      | {z} = ˜ B w w k = ˜ A ˜ x k + ˜ B w w k (7) for k = 0 , 1 , . . . , which is an L TI model that we will refer to as the p - augmented L TI system. It should also be noted that, if we consider a time-varying FOS system ∆ α x k +1 = A k x k + w k , (8) then (3) would still be valid, except that A j would need to be index ed by k as well and we could use the same reasoning to deri ve a finite-history approximation ˜ x k +1 = ˜ A k ˜ x k + ˜ B w w k , (9) with the only exception that this p -augmented model is linear and time-varying . These representations will play a key role in the imple- mentation of our proposed model predictiv e control strat- egy , since such an MV AR approximation of an FOS plant, re-written as an L TI model, allows for ef ficient numerical solutions to be determined for quadratic-cost optimal control. The reason for this is that a full representation of an FOS model would require increasing memory storage, while making the problem of computing optimal control actions an increasingly computationally demanding task. Naturally , such a completely faithful representation of fractional-order systems would therefore become intractable for most practi- cal applications at likely little cost benefit, since the temporal dependence, while a long-term one for FOS models, is also decaying with respect to the weights. A. Model Pr edictive Contr ol W e no w focus on the design of a full-state feedback controller for a linear time-v arying system ov er discrete-time, x k +1 = A k x k + B k u k + B w k w k , (10) where w k denotes a sequence of independent and identi- cally distributed (i.i.d.) random vectors, following a standard normal distribution with zero mean and identity cov ariance matrix. The weight matrices B w k are intended to make the notation cleaner and add flexibility to the process noise w 0 k = B w k w k since the cov ariance matrix of w 0 k does not need to be positi ve definite, which is a strict requirement for Gaussian distributions. The objective is to design the feedback controller such that it minimizes a quadratic cost functional of the input and state vectors over an infinite time horizon. In other words, the objecti ve is to determine the sequence of control inputs u 0 , . . . , u N − 1 that minimizes a quadratic cost function of the form minimize u 0 ,...,u N − 1 E ( N X k =1 k x k k 2 Q k + N X k =1 c T k x k + N − 1 X k =0 k u k k 2 R k ) subject to x 0 = x init , x k +1 = A k x k + B k u k + B w k w k , for k = 0 , 1 , . . . , N − 1 , (11) for N → ∞ , where x k ∈ R n , u k ∈ R n u , and Q 1 , . . . , Q N ∈ R n × n and R 0 , . . . , R N − 1 ∈ R n u × n u are giv en positi ve semidefinite matrices. Recall that, Q ∈ R n × n is a positive semidefinite matrix if x T Qx ≥ 0 for every x ∈ R n , and k x k Q = p x T Qx in that case. For seizure mitigation via TMS, we propose to use Q k = I n × n , c = 0 n × 1 , and R k = εI n u × n u with ε > 0 , such that the objectiv e becomes largely to steer the total energy in the expected value of the brain signals tow ards the smallest amount possible. The quadratic term on the input, which represents the stimulation signal induced in the brain via TMS, is intended to add a penalization term for stimulating the patient too harshly , since this may be unsafe, create discomfort for the patient, or hav e harmful psychological effects [35]. Problem (11) is often referred to as the linear-quadr atic-Gaussian (LQG) control problem [36], and it may be solved offline and implemented in real time through a feedback control strategy of the form u k = K k x k ( k = 0 , 1 , . . . , N − 1) , (12) for some appropriate matrices K 0 , . . . , K N − 1 ∈ R n u × n (called feedback gain matrices) [37], which are functions of x 0 . Furthermore, (11) can be written as an unconstrained quadratic pr ogram (QP), i.e., the minimization of a function f ( x ) = 1 2 k x k 2 Q + c T x . Such problems can be efficiently solved numerically , for instance, via the quadprog() func- tion in M A T L A B . In fact, even if we include additional linear constraints on the state and input vectors, the problem will remain a QP . In particular , we may consider constraints of the form u min ≤ ( u k ) i ≤ u max for e very time step k and ev ery input component indexed by i . This can be used to efficiently ensure that our proposed TMS approach will only ev er administer safe v oltage stimuli. In model predictiv e control, the objectiv e is to solve these consecutiv e constrained finite-horizon LQG problems over a moving horizon, in order to create additional robustness, as compared to directly solving for N → ∞ . More precisely , at time step k , the proposed strategy inv olves solving the problem minimize u k ,...,u k + P − 1 E ( P X j =1 k x k + j k 2 Q k + j + P X j =1 c T k + j x k + j + P − 1 X j =0 k u k + j k 2 R k + j ) subject to x k = observed or estimated current state , x k + j +1 = A k + j x k + j + B k + j u k + j + B w k + j w k + j , for j = 0 , 1 , . . . , P − 1 , other linear constraints on x k +1 , . . . , x k + P , u k , . . . , u k + P − 1 , (13) where P is called the pr ediction horizon , but it only deploys the control strategy associated with the first M time steps (referred to as the contr ol horizon ). Simply speaking, after we reach the state x k + M − 1 , we update k ← k + M − 1 and recompute the new LQG solution. This way , by design, there is a new layer of rob ustness that the solution of an LQG (ev en for N → ∞ ) does not offer by itself, since the optimal strategy is constantly being re-ev aluated based on short-term control action implementation of a long-term prediction [38]. I V . S I M U L A T I O N R E S U LT S In what follows, we propose to illustrate the use of the fractional-order model predicti ve control frame work for TMS in the context of annihilating or mitigating epileptic seizures using real EEG data. W e start by considering an epileptic seizure captured by a linear FOS model whose parameters are obtained through a system identification method using real brainwave data. Then, we simulate three dif ferent stim- ulation strategies: (i) an open-loop stimulation strategy; (ii) an ev ent-triggered open-loop strategy (i.e., once a certain ev ent, in this case the beginning of a seizure, is detected, a short-term open-loop stimulation strategy is deployed); and (iii) a closed-loop stimulation strategy based on our proposed FOS-MPC approach. System Identification Thr ough P arameter Estimation with EEG data First, we need to determine the parameters A and α that model both spatial coupling and fractional exponents, respectiv ely , that craft the e volution of the state x k ∈ R n in the FOS model ∆ α x k +1 = Ax k + w k , (14) with w k denoting additiv e white Gaussian noise (A WGN). W e will consider only n = 4 components in our state v ector x k with each one corresponding to the measurement obtained from a single EEG microelectrode. T o identify the parameters A , α , and the cov ariance matrix Σ w (assumed equal to the diagonal matrix Σ = σ 2 w I 4 × 4 for simplicity) associated with w k , we used the method proposed in [16]. For illustration purposes, we consider the normalized (i.e., voltage units ha ve been scaled to the interv al − 1 to 1 ) EEG recordings 1-4 of subject 11 from the CHB-MIT Scalp EEG database [39]. The parameters obtained are as follows: A =     0 . 0402 0 . 0604 − 0 . 0040 − 0 . 0450 0 . 0340 − 0 . 0571 0 . 0742 0 . 0701 − 0 . 0119 − 0 . 0032 − 0 . 0105 0 . 0078 − 0 . 0335 0 . 0165 − 0 . 0009 0 . 0453     , (15a) α =  0 . 6606 0 . 7973 1 . 0670 0 . 6977  T , (15b) and σ 2 w = 0 . 2 . More precisely , we utilized a normalized 10-second sample during a period of ictal acti vity (i.e., activity during an epileptic seizure), sampled at 160 Hz. Experiment 1: Open-Loop Electrical Neur ostimulation Based upon the identified parameters of the system (14)–(15), we consider the follo wing forced FOS ∆ α x k +1 = Ax k + B u k + w k + d k , (16) where d k denotes a disturbance term that can be understood as a persistent neural activity incoming from nearby regions. For the sake of our simulations, d k will be set as frequent wa velet-lik e bursts of amplitude d k = 0 . 25 and d k = 1 , which occur in a disjoint manner at random points of time according to Poisson counting processes with rates of 0.2 s and 1 s, respectively . Fig. 1. Simulated seizure is depicted (in blue) versus an attempted controlled signal (in red) under an open-loop stimulation strategy (in yellow). Additionally , u k denotes the voltage stimulus signal be- ing applied. For the sak e of simplicity , we will assume B = [1 , 1 , 1 , 1] T , which corresponds to a stimulus that per- turbs all the channels uniformly . Sev eral open-loop stimulation strategies were imple- mented, which always entail u k to be determined without information about the state of the system. Specifically , when considering the previous model, we start actuation at the 4 second mark, when the EEG signal (blue) branches out into two signals, one being the continuation of the original signal without an y stimulation, and the other being the electrically stimulated version (red) of the signal through the open-loop controlled stimulus (orange). Here, the input will consist of 1 s periods of sinusoidal activity of amplitude 0.5, frequency 16 Hz, phase 0, and consecutiv ely followed by unstimulated periods lasting 0.5 s. As we can see in Fig. 1, this strategy is unable to steer the ev olution of the system tow ards normal activity . Specifically , the controlled brainwav es seem to be almost unaf fected, which is consistently observed if we adopt slightly different open-loop strategies (e.g., adopting a different wav eform, as it is in the case of a biphasic rectangular pulse). In fact, we have observed that such open-loop control strategies can induce the increase of activity and, at times, even originate seizure-like activity . As a particular instance of these findings, we re-identified the parameters in our system based on inter-ictal data of the same subject as before. Next, we considered a stimulation strategy as depicted in Fig. 2, which lead to sev eral periods of seizure-like activity which would not hav e occurred if the original system was left to ev olve on its own. Fig. 2. Seizure-like activity (in red) is induced by an open-loop stimulation strategy (in yellow) of an otherwise regular activity depicted in blue. Experiment 2: Event-T rigger ed Open-Loop Contr ol Whereas in the Experiment 1 we focused on open-loop strategies that start at fully pre-specified instances of time, in practice, neurostimulation devices such as those used in TMS implement the so-called ev ent-triggered open-loop stimulation strategy . More specifically , the (open-loop) stim- ulation strategy is only activ ated when a phenomenon of interest is observed that can be described by the previously obtained data. In what follo ws, we assumed that, when activ ated due to the detection of seizure-like acti vity , the same open-loop stimulation strategy as presented in Fig. 2 with twice the amplitude is deployed – see Fig. 3. In practice, there are a variety of such seizure-like activity detectors that can be considered (e.g., the line length feature that measures the length of the line described by the activity during a windo w of time [40]). Notwithstanding, we obtained a similar conclusion as when open-loop strategies were not triggered by an event. Fig. 3. Simulated seizure is depicted (in blue) versus an attempted controlled signal (in red) under an event-triggered open-loop stimulation strategy (in yellow). Experiment 3: Closed-Loop Electrical Neur ostimulation us- ing MPC on FOS Finally , we implemented our proposed TMS strategy using MPC on the FOS models described earlier in this paper . More precisely , for the cost function in (13), we utilized Q k = 10 I n , R k = I n u , and c k = 0 n u × 1 , to emphasize mini- mizing the overall energy in the measured brainwav es, while penalizing slightly for overly aggressi ve stimulation. Further- more, we included a safety linear constraint of − 1 ≤ u k ≤ 1 . Our predictiv e model was based on a ( p = 4) − step (2.5 ms) MV AR predictiv e model approximation of the FOS plant, with a ( P = 32) − step (20 ms) prediction horizon and ( M = 8) − step (5 ms) control horizon. Fig. 4. Simulated seizure is depicted (in blue) versus controlled signals (in red) for all channels, under our proposed closed-loop stimulation strategy (in yellow). The results presented in Fig. 4 include all 4 controlled channels controlled simultaneously , since failure to driv e any one of them tow ards a normal range would imply failure in the seizure mitigation objecti ve. W e can see that the proposed strategy using TMS achie ves the desired goal and implicitly provides us with a detector , giv en that it tends to provide virtually no stimulation except during very brief periods of time, at which point only (time-varying) impulse-like stimuli are deployed. Remark 1. Any control scheme that requires computa- tionally demanding large-scale optimization methods to be in volved, will, in turn, require some form of approximation in order to increase computational efficiency and to enable true real-time control of the system. As such, it is imperativ e to understand the trade-offs in computational performance in volved when using such approximations, which naturally depend from system to system. Remark 2. The actual dynamics of the brain are highly nonlinear and time-varying. For this reason, it is crucial to re-identify the fractional-order parameters in our predictiv e model in a real-life implementation of our proposed strategy . An empirical analysis of parameter tuning via trial and error suggests that an L TI predictiv e model approximation of the FOS plant with non-augmented state (i.e., an MV AR( p ) predictiv e model with p = 1 ) is inadequate to achie ve satisfactory lev els of performance in the context of seizure mitigation. Howe ver , it also suggests that the memory of the predictiv e model can be limited to p = 8 past time steps (ranging 0.5 ms) for the conducted experiments. Anything beyond that led to negligible gains in performance at a high computational cost. Remark 3. It should also be noted that in the proposed FOS-MPC stimulation strategy , there are still some design parameters that need to be manually calibrated, such as the prediction horizon P , the control horizon M , the memory horizon p , and the input energy penalization weight ε ≥ 0 . Notwithstanding these considerations, there is a considerable theoretical foundation dedicated to studying the design of MPC algorithms that achiev e stability , rob ustness, and other performance guarantees [41]–[44]. This body of results may be used to guide and systematize the parameter calibration stage under a sound and justifiable basis. V . C O N C L U S I O N S A N D F U T U R E W O R K W e presented a methodological framework towards real-time feedback control with constraints for neurophys- iological systems. Specifically , we pedagogically introduced a model predicti ve control (MPC) approach when the neuro- physiological process can be modeled by a fractional-order system. In doing so, we focused on TMS for epilepsy , and using systems with seizure-like characteristics, we showed that the stimulation strategies obtained by the proposed framew ork enabled us to annihilate and mitigate epileptic seizures. Although we ha ve focused mainly on TMS for epilepsy , we believ e that the proposed frame work can be readily applied to other forms of neurostimulation with an adequate change in the optimization problem (i.e., in the objectiv e function and constraints). R E F E R E N C E S [1] N. Brodu, F . Lotte, and A. L ´ ecuyer , “Exploring two novel features for EEG-based brain–computer interfaces: Multifractal cumulants and predictiv e complexity , ” Neur ocomputing , vol. 79, pp. 87–94, 2012. [2] E. A. Ihlen and B. V ereijken, “Interaction-dominant dynamics in hu- man cognition: Beyond 1 /f ( α ) fluctuation. ” Journal of Experimental Psychology: General , vol. 139, no. 3, p. 436, 2010. [3] P . Ciuciu, G. V aroquaux, P . Abry , S. Sadaghiani, and A. Kleinschmidt, “Scale-free and multifractal properties of fMRI signals during rest and task, ” Fr ontiers in Physiology , vol. 3, p. 186, 2012. [4] T . Zorick and M. A. Mandelkern, “Multifractal detrended fluctuation analysis of human EEG: preliminary investigation and comparison with the wavelet transform modulus maxima technique, ” PloS one , vol. 8, no. 7, p. e68360, 2013. [5] Y . Zhang, W . Zhou, and S. Y uan, “Multifractal analysis and rele- vance vector machine-based automatic seizure detection in intracranial EEG, ” International Journal of Neural Systems , v ol. 25, no. 06, p. 1550020, 2015. [6] D. Papo, “Functional significance of complex fluctuations in brain activity: from resting state to cogniti ve neuroscience, ” F r ontiers in Systems Neuroscience , vol. 8, p. 112, 2014. [7] J. Suckling, A. M. Wink, F . A. Bernard, A. Barnes, and E. Bull- more, “Endogenous multifractal brain dynamics are modulated by age, cholinergic blockade and cognitiv e performance, ” Journal of Neur oscience Methods , vol. 174, no. 2, pp. 292–300, 2008. [8] V . E. T arasov , Fr actional Dynamics: Applications of F ractional Cal- culus to Dynamics of P articles, Fields and Media . Springer Science & Business Media, 2011. [9] Y . Xue, S. Rodriguez, and P . Bogdan, “ A spatio-temporal fractal model for a CPS approach to brain-machine-body interfaces, ” in Pr oceedings of the Design, Automation & T est in Europe Confer ence & Exhibition (D A TE), 2016 . IEEE, 2016, pp. 642–647. [10] Y . Xue, S. Pequito, J. R. Coelho, P . Bogdan, and G. J. Pappas, “Minimum number of sensors to ensure observability of physiological systems: A case study , ” in Pr oceedings of the 54th Annual Allerton Confer ence on Communication, Control, and Computing . IEEE, 2016, pp. 1181–1188. [11] Y . Xue and P . Bogdan, “Reliable multi-fractal characterization of weighted complex networks: algorithms and implications, ” Scientific Reports , vol. 7, no. 1, p. 7487, 2017. [12] ——, “Constructing compact causal mathematical models for complex dynamics, ” in Proceedings of the 8th International Confer ence on Cyber-Physical Systems . A CM, 2017, pp. 97–107. [13] S. Pequito, P . Bogdan, and G. J. Pappas, “Minimum number of probes for brain dynamics observability , ” Proceedings of the 54th IEEE Conference on Decision and Contr ol (CDC) , pp. 306–311, Dec 2015. [14] S. Pequito, A. Clark, and G. J. Pappas, “Discrete-time fractional-order multiple scenario-based sensor selection, ” in Proceedings of the 2017 American Control Conference (A CC) , May 2017, pp. 5488–5493. [15] V . Tzoumas, Y . Xue, S. Pequito, P . Bogdan, and G. J. Pappas, “Select- ing sensors in biological fractional-order systems, ” IEEE T ransactions on Control of Network Systems , 2018. [16] G. Gupta, S. Pequito, and P . Bogdan, “Dealing with Unknown Un- knowns: Identification and Selection of Minimal Sensing for Frac- tional Dynamics with Unknown Inputs, ” in Proceedings of the 2018 American Control Conference (A CC) , June 2018, pp. 2814–2820. [17] ——, “Re-thinking EEG-based non-inv asive brain interfaces: model- ing and analysis, ” in Proceedings of the 9th ACM/IEEE International Confer ence on Cyber-Physical Systems . IEEE Press, 2018, pp. 275– 286. [18] A. Dzielinski and D. Sierociuk, “ Adaptive feedback control of frac- tional order discrete state-space systems, ” in CIMCA-IA WTIC , 2005. [19] E. F . Camacho and C. B. Alba, Model Predictive Contr ol . Springer Science & Business Media, 2013. [20] D. Q. Mayne, J. B. Rawlings, C. V . Rao, and P . O. Scokaert, “Constrained model predictiv e control: Stability and optimality , ” Au- tomatica , vol. 36, no. 6, pp. 789–814, 2000. [21] S. J. Qin and T . A. Badgwell, “ A survey of industrial model predictive control technology , ” Control Engineering Practice , vol. 11, no. 7, pp. 733–764, 2003. [22] F . Allg ¨ ower and A. Zheng, Nonlinear Model Predictive Contr ol . Birkh ¨ auser , 2012, vol. 26. [23] C. E. Garcia, D. M. Prett, and M. Morari, “Model predictive control: theory and practice–a surve y, ” Automatica , v ol. 25, no. 3, pp. 335–348, 1989. [24] M. Morari and J. H. Lee, “Model predictive control: past, present and future, ” Computers & Chemical Engineering , vol. 23, no. 4-5, pp. 667–682, 1999. [25] W . M. Schuepbach, J. Rau, K. Knudsen, J. V olkmann, P . Krack, L. Timmermann, T . H ¨ albig, H. Hesekamp, S. Nav arro, N. Meier , et al. , “Neurostimulation for Parkinson’ s disease with early motor complications, ” New England Journal of Medicine , vol. 368, no. 7, pp. 610–622, 2013. [26] K. Witt, O. Granert, C. Daniels, J. V olkmann, D. Falk, T . van Eimeren, and G. Deuschl, “Relation of lead trajectory and electrode position to neuropsychological outcomes of subthalamic neurostimulation in Parkinson’ s disease: results from a randomized trial, ” Brain , vol. 136, no. 7, pp. 2109–2119, 2013. [27] R. Nardone, Y . H ¨ oller , F . T ezzon, M. Christov a, K. Schwenker, S. Go- laszewski, E. Trinka, and F . Brigo, “Neurostimulation in Alzheimer’ s disease: from basic research to clinical applications, ” Neurolo gical Sciences , vol. 36, no. 5, pp. 689–700, 2015. [28] L. Marangell, M. Martinez, R. Jurdi, and H. Zboyan, “Neurostim- ulation therapies in depression: a re view of new modalities, ” Acta Psychiatrica Scandinavica , vol. 116, no. 3, pp. 174–181, 2007. [29] B. H. Bewernick, R. Hurlemann, A. Matusch, S. Kayser, C. Grubert, B. Hadrysiewicz, N. Axmacher , M. Lemke, D. Cooper-Mahkorn, M. X. Cohen, et al. , “Nucleus accumbens deep brain stimulation decreases ratings of depression and anxiety in treatment-resistant depression, ” Biological Psychiatry , vol. 67, no. 2, pp. 110–116, 2010. [30] V . Sturm, D. Lenartz, A. Koulousakis, H. T reuer, K. Herholz, J. C. Klein, and J. Klosterk ¨ otter , “The nucleus accumbens: a target for deep- brain stimulation in obsessive-compulsi ve and anxiety disorders, ” in Pr oceedings of the Medtronic F orum for Neur oscience and Neuro- T echnology 2005 . Springer, 2007, pp. 62–67. [31] J.-P . Lefaucheur , N. Andr ´ e-Obadia, A. Antal, S. S. A yache, C. Baeken, D. H. Benninger, R. M. Cantello, M. Cincotta, M. de Carvalho, D. De Ridder , et al. , “Evidence-based guidelines on the therapeutic use of repetitive transcranial magnetic stimulation (rTMS), ” Clinical Neur ophysiology , vol. 125, no. 11, pp. 2150–2206, 2014. [32] J. Fang, M. Zhou, M. Y ang, C. Zhu, and L. He, “Repetitive transcranial magnetic stimulation for the treatment of amyotrophic lateral sclerosis or motor neuron disease, ” Cochrane Database of Systematic Reviews , no. 5, 2013. [33] M. A. Dimyan and L. G. Cohen, “Contribution of transcranial mag- netic stimulation to the understanding of functional recovery mecha- nisms after stroke, ” Neur orehabilitation and Neural Repair , v ol. 24, no. 2, pp. 125–135, 2010. [34] S. Guermah, S. Djennoune, and M. Bettayeb, “Controllability and observability of linear discrete-time fractional-order systems, ” Interna- tional Journal of Applied Mathematics and Computer Science , vol. 18, no. 2, pp. 213–222, 2008. [35] S. Moratti and D. Patterson, “ Adverse psychological ef fects to deep brain stimulation: Overturning the question, ” AJOB Neur oscience , vol. 5, no. 4, pp. 62–64, 2014. [36] H. Kwakernaak, Linear Optimal Control Systems , R. Sivan, Ed. New Y ork, NY , USA: John W iley & Sons, Inc., 1972. [37] J. P . Hespanha, Linear systems theory . Princeton University Press, 2018. [38] B. W . Bequette, “ Algorithms for a closed-loop artificial pancreas: The case for model predicti ve control, ” Journal of Diabetes Science and T echnology , vol. 7, no. 6, pp. 1632–1643, 2013. [39] A. L. Goldberger , L. A. N. Amaral, L. Glass, J. M. Hausdorff, P . C. Ivano v , R. G. Mark, J. E. Mietus, G. B. Moody , C.-K. Peng, and H. E. Stanley , “PhysioBank, PhysioT oolkit, and PhysioNet: Components of a new research resource for complex physiologic signals, ” Circulation , vol. 101, no. 23, pp. e215–e220, 2000 (June 13). [40] R. Esteller, J. Echauz, T . Tcheng, B. Litt, and B. Pless, “Line length: an efficient feature for seizure onset detection, ” in Pr oceedings of the 23r d Annual International Confer ence of the IEEE Engineering in Medicine and Biology Society , vol. 2. IEEE, 2001, pp. 1707–1710. [41] X. Cheng and B. H. Krogh, “Stability-constrained model predictiv e control, ” IEEE T ransactions on A utomatic Contr ol , vol. 46, no. 11, pp. 1816–1820, Nov 2001. [42] N. H. El-Farra, P . Mhaskar, and P . D. Christofides, “Uniting bounded control and MPC for stabilization of constrained linear systems, ” in Pr oceedings of the 2002 American Control Conference , vol. 6, May 2002, pp. 4493–4500 vol.6. [43] M. N. Zeilinger , C. N. Jones, D. M. Raimondo, and M. Morari, “Real- time MPC - stability through robust MPC design, ” in Proceedings of the 48h IEEE Conference on Decision and Contr ol (CDC) held jointly with 2009 28th Chinese Control Confer ence , Dec 2009, pp. 3980–3986. [44] J. Pannek and K. W orthmann, “Stability and performance guarantees for model predictive control algorithms without terminal constraints, ” ZAMM Journal of applied mathematics and mechanics: Zeitschrift fr angewandte Mathematik und Mechanik , v ol. 94, pp. 317–330, 04 2014.