Modeling and Control of Biological Systems: Applications in Tissue Engineering and Artificial Pancreas

Abstract

The complexity of biological systems demands the use of appropriate models in order to relate the effects of inputs on the system and capture events that influence the final quality of final products. Several modeling paradigms can be used to describe system behavior. Agent-based modeling provides an efficient and powerful environment to represent characteristics at organ, cellular, and subcellular levels, enables interactions of various components of a system based on events triggered during a simulation, and tracks dynamic information of individual entities. Agents are individual building blocks that are capable of performing various actions independently based on both their own perception of their environment, and the quantitative and qualitative rules that govern their behavior. Three-dimensional tissue engineering where angiogenesis and tissue growth occur in a scaffold illustrates agent-based modeling and its use for supervision of tissue growth, prediction of the final tissue properties, and regulation of the growth process for the desired final tissue. Time series models offer an alternative modeling paradigm and generate input-output (blackbox) models that have been successful in supervision and control of complex processes. A multivariable control system for an artificial pancreas based on recursive time series models and generalized predictive control illustrates the performance of this control system for regulating glucose concentration in patients with type 1 diabetes.

Introduction

The potential impact of simulation-based engineering science on advances in science and technology and current challenges are documented in a recent NSF report.¹ The challenges include difficulties inherent in multi-scale modeling, the design and implementation of dynamic data-driven application systems, the improvement of methods to quantify uncertainty and model validation, and better integration of data-intensive computing, visualization, and simulation.¹ Expected payoffs are dramatic advances on a broad front including medicine and materials science.¹ Many biological systems are the result of interactions of simpler behaviors, yielding a complex system that defies understanding through intuition and “it is critical to have a model that helps us understand the structure of these phenomena.”² In contrast with these models that describe physical and biological phenomena, input-output (black box) models are more useful in developing model-based supervision and control systems. They have been used successfully in pharmaceutical manufacturing and medical devices. This paper presents modeling techniques for both types of activities, agent-based modeling for describing three-dimensional (3D) tissue engineering where angiogenesis and tissue growth occur in a scaffold, and recursive time series modeling for model-based control of artificial pancreas systems.

Agent-Based Modeling

Agent-based modeling (ABM) provides an efficient and powerful environment for representing the characteristics and describing the behavior of a system and its components at various levels and tracking dynamic information of individual entities. Agents are individual building blocks that are capable of performing various actions proactively and independently based on their perception of their environment, and the quantitative and qualitative rules that govern their behavior. A critical feature of ABM is the instantaneous response of agents to changes taking place in their environment during a simulation at each time step of execution. An agent senses the changes in its environment and neighboring agents (cells in our model) as they occur, and makes instantaneous decisions about its own behavior. The ABM approach can track system behavior at various levels including changes in individual cells and also provide population statistics.

An agent is an intelligent, autonomous, and interacting software entity that represents an active constituent of a system. It makes decisions based on the analysis of the conditions and expectations in its environment and its internal logic, and performs actions based on these decisions. The behavior of the system emerges from the simple, lower-level actions of its constituents. Software agents utilize their interaction and collaboration abilities to perform actions while their environments change in predictable or unpredictable manners. The characteristics of agents include³:

Autonomy: The ability of an agent to make decisions based on its rules and characteristics, and environmental conditions. For example, one cell agent can survive, grow, migrate, or proliferate in a changing environment while its neighbor agent may perform a different action based on stochastic variations in its attributes and local conditions.

Communication: An agent can communicate with other agents and its environment. For example, a cell agent can receive information about metabolite, nutrient, and oxygen concentrations, discern the presence of cells and blood vessels in its vicinity, and send signals such as growth factors (GF) about its own status.

Cooperation: An agent can collaborate and cooperate with other agents in order to achieve a common goal. For example, agents representing endothelial cells (EC) form blood vessels and anastomosis (the connection of 2 vessel branches leading to blood perfusion).

Agents can also cooperate to solve a problem such as growing blood vessels in the direction of tissue cells in need of nutrients and oxygen.

Mobility: An agent can physically move in its environment at its own discretion. For example, tissue cell agents migrate in the direction of higher concentrations of nutrients and oxygen (in the direction of functional vasculature). EC agents can also become fixed at a location. For example, the locations of stalk cells of established blood vessels are fixed.

ABM becomes more appealing when the system to be modeled becomes increasingly complex and heterogeneous, providing model features that conventional equation-based methodologies cannot describe efficiently. ABMs can describe and forecast the dynamics of system components by tracking the behavior of each agent and enabling the emergence of higher-scale properties from the rules that govern the lower-scale system components. New hypotheses can be tested by changing the rules, relations and attributes of agents, leading to the identification of experimental conditions to investigate promising hypotheses that may lead to desirable phenomena and product properties. ABMs can be maintained and refined easily because of their modular and object-oriented structure by changing rules and attributes, and comparing simulation results with observed phenomena.

Biological systems are multi-scale and heterogeneous complex systems, with many types of cells and metabolites, and spanning different scales from molecular to cellular and tissue level. ABMs are suitable for representing such complex biological systems. The agents can represent organs, tissues, cells, or genes. ABMs have been used successfully in recent years for modeling biological systems focusing on individual cell behavior such as angiogenesis,^4,5 osteogenesis of stem cells,⁶ brain tumors,⁷ and arterial adaptation in hypertension.⁸ In another work, an ABM coupled with a transport model to simulate the dynamic relationship between a cell and its microenvironment was developed for a tissue engineering bioreactor.⁹

Agent-based modeling of angiogenesis and tissue growth

We developed a multi-layer ABM to simulate the process of vascularized bone regeneration of multi-potent stem cells in 3D porous polymeric scaffolds from stem cells and chemical factors to optimize functional tissue formation. The ABM software consists of 3 interacting layers (Figure. 1). One layer describes angiogenesis, focusing on the behavior of individual ECs to form vessel sections in response to GF gradient in the environment. EC agents mimic the actions that ECs can make during angiogenesis, including elongation, proliferation, tip cell migration, and anastomosis, as governed by a rule base.^4,10 The second layer represents the scaffold structure and describes events occurring within the 3D space, including the location of the scaffold, diffusion of soluble factors to promote angiogenesis and bone regeneration. Dynamic gradients of vascular endothelial growth factor (VEGF) and bone morphogenic protein-2 (BMP-2; for the formation of bone) release data are incorporated into the model in the scaffold layer. The third layer describes the behavior of multi-potent stem cells (MSCs) with osteogenic potential and governs the differentiation to osteoblast (OB) and the formation of bone tissue. The main agent of this layer is the stem cell agent.⁶

Figure 1. The interacting layers in the multi-layer ABM. Tissue growth layer includes MSC and OB agents; angiogenesis layer includes EC agent and scaffold-soluble layer includes the porous scaffold and the release concentration data of VEGF and BMP-2.

Agents and Rule-base

A set of rules developed by harvesting the literature is collected in a rule base to govern the behavior of individual EC and MSC agents. The rule base constitutes the internal logic of agents and controls the sequential actions performed by each cell during simulation. The states and the actions that agents can take are summarized in Table 1.

Table 1. Actions and states of agents in angiogenesis and tissue growth system

EC agent rule-base leads to actions including migration, elongation, proliferation, branching, and connection of blood vessels (anastomosis). An activated tip cell initially schedules the occurrence of the next migration action. The time interval between migration actions depends upon the ligand density in the agent’s environment. Next, the tip cell checks the “persistency time,” and if it is passed, migration direction is updated as the direction of highest GF concentration gradient. After completing the migration, the tip cell searches its environment, and if another EC agent is in its vicinity, anastomosis takes place. Stalk EC elongation is permitted until it reaches “maximum EC length.” Once the length of an EC reaches its maximum, proliferation occurs and another EC is generated. The newly generated EC agent becomes the leading node (tip cell) and starts migration while the parent agent becomes fixed (stalk cell), and is no longer able to migrate further.

Another set of rules is derived for the tissue layer and its agents (MSCs, OBs). An MSC agent is capable of sensing the growth factor concentration in its vicinity, migrating towards the BMP-2 gradient, going through cell cycle, growing in size, contributing to the VEGF concentration in the environment and committing and then differentiating to a bone cell (OB agent) in a continuous space. Two levels are built in the rule set that govern the behavior of MSC agents. The first level is the base level rule set to regulate the cell cycle and survival of cells. A higher-level rule set is introduced to the model to update the parameters of cells going through the differentiation. In the cycle rule set, the first thing each agent does is to calculate the distance to the closest stable blood vessel. Only the stable (anastomosed) blood vessels are considered as having blood flow and providing a source for oxygen and nutrition. If there is no blood vessel in close vicinity (0 to 100 μm), the cell will change its state to quiescent phase (G₀). Hypoxic cell state (H) will be assigned to the cell if there is no blood vessel in a further distance (100 to 200 μm) and start the hypoxia clock. If no blood vessel reaches the close vicinity of a hypoxic cell, the cell will go through apoptosis. If the cell is in the vicinity of a stable blood vessel, then it will go through the cell cycle and spend the required time at each state and complete the actions before it could change the state in the cycle. Cells can only proliferate at the end of the G₂M phase. Both MSCs and OBs migration are governed by chemotactic random migration through BMP-2 gradient. Cells can differentiate only in G1 state, and this is when the higher-level rule set is called. The “Differentiation” rule set is responsible for changing the parameters in the “Cycle” rule set based on the environmental condition (BMP-2 concentration). Three snapshots of a simulation illustrate predicted angiogenesis and tissue growth resulting from agent actions over time in Figure 2.

Figure 2. Angiogenesis and tissue growth in cell-seeded porous scaffolds (white). Different colors represent different phases in cell cycle.

Time Series Modeling

Time series (TS) models have been successfully used for modeling the effects of inputs on outputs in many areas ranging from economic forecasts to process operations. Linear TS models can be represented as:

where k denotes the sampling time, y(k) is the system output, u_i(k) is the ith input variable, Єk is the white noise, di is the delay term for the corresponding input. The polynomials that represent the past values of the inputs and outputs are expressed by using the backward shift operator q^-1

where n_A, n_Bi, and n_C are model orders. The model Eq (1) is called multivariable autoregressive moving average model with exogenous inputs (ARMAX) since it includes the past values of outputs (autoregressive terms), the past values of inputs adjusted for appropriate delays (exogenous terms) and a random noise sequence that accounts for stochastic variations (moving average terms) to compute the current value y(k).

A single linear model cannot describe all the states and dynamic variations of complex systems such as glucose concentration (GC) dynamics in humans in response to meals, physical activity, or stress. Yet the simple structure of the TS model is very appealing for use in a model-based control system. A recursive TS model that updates the model’s parameters with every new set of measurements provides an effective solution. Given a time series model, weighted recursive least squares (WRLS) method¹¹ estimates the model parameters by minimizing the error between the observed and estimated values of the output(s) in the sliding time window.

Writing the ARMAX model in linear regression form:

The white noise term in Eq (6) is replaced with model prediction error (residual) e(k) when the model is used for predictions since the former is unknown.¹²

Optimization methods may yield unstable recursive models for a stable system because of measurement noise or use of open-loop model identification, estimation, and diagnostic criteria with closed-loop data. We developed a constrained optimization method that uses the state space model extracted from the TS model for determining the parameters of the recursive model that assures model stability.¹³

The unknown model parameters in Eq (7) are calculated by solving of a constrained optimization problem:

where P(k) is the estimate of the error covariance matrix and λ is the forgetting factor that adjusts the weight of recent measurements with respect to older ones. This enables faster convergence to new parameter values when significant changes occur in a patient’s glucose level. The model parameters are updated recursively after receiving new measurements at each sampling time. The state space model is asymptotically stable if all eigenvalues of state matrix are inside the unit circle. One way to ensure the stability constraints is to impose the condition that the spectral radius ρ(Ã) of the state matrix à be <1.¹³

Multivariable Artificial Pancreas Control System

Patients with type 1 diabetes need to administer exogenous insulin to maintain their blood glucose levels in a target range. An artificial pancreas (AP) system calculates the appropriate amount of insulin to be infused by an insulin pump based on glucose concentration measurements from a continuous glucose monitoring (CGM) system. Significant advances have been made in the development of an AP in recent years.^12-20 Glycemic control can be achieved by a completely automated external closed-loop insulin delivery. Various automatic control strategies with different modeling algorithms have been investigated during the last 3 decades.^14-19 Most of the AP systems are based on CGM-only measurements and users often manually enter meal information and other important changes such as exercise or sleep to improve controller performance.^14-19 However, additional biometric/ physiological measurements can provide valuable information that can eliminate the need for manual entries. An adaptive control system that relies on subcutaneous glucose concentration measurements from a CGM and biometric/physiological data from a sports armband provides a powerful alternative that does not necessitate any meal or activity announcements.²⁰ We developed an adaptive control system based on the generalized predictive control (GPC) approach by using multivariable recursive time series models. Model parameters are updated with every new measurement set collected at each sampling time (which was 10 minutes in our system).

Physiological signals such as energy expenditure and galvanic skin response are used along with glucose measurements (collected via CGM) to generate a multiple-input-single-output model for predicting future glucose concentrations used by the controller.

After identification of unknown parameters in Eq (7), the ARMAX model is used in the GPC for calculation of the appropriate insulin amount to be infused by a constrained minimization problem:

where r(∙) is the target glucose trajectory and w(∙) is the adaptively calculated weight matrix.¹² u_min and u_max are the minimum and maximum limits for calculated insulin and N₁ and N₂ are the first and last time instants of the prediction horizon and N_u is the control horizon. In addition to using physiological signals as inputs of the model in Eq (1) for prediction of glucose concentration, the current values of physiological signals are used in a penalty function of the controller to reduce the insulin infusion rate during exercise. The insulin-on-board is also estimated and used in control decisions. The control system is integrated with a hypoglycemia early alarm system that predicts potential hypoglycemia episodes 30 minutes in advance and warns the patient to take a rescue carbohydrate (Figure 3).¹³

Figure 3. Integrated multivariable adaptive artificial pancreas control and hypoglycemia early alarm system.

The integrated multivariable control and hypoglycemia prediction alarm system was tested in clinical studies that lasted for 60 hours without any meal or activity announcements.¹⁹ The adaptive control system kept glucose concentration in the normal pre-prandial and post-prandial range (70 to 180 mg/dL) during most of the test period without any meal or activity announcements (Figure 4). The subjects performed medium-intensity exercise every day. They started by walking on a treadmill using ramped protocol to transition to a run where the speed (miles per hour [mph]) and the incline (%) of the treadmill were gradually increased until the subjects terminated the exercise session based on self-reported symptoms (eg, exhaustion). Overall, the subjects exercised 85 ± 8.3% (Mean and Standard Deviation) of their age-predicted heart rate of 220-age (range: 70% to 97%). Glucose concentration (CGM values and infrequent blood glucose measurements), insulin infusion rates and energy expenditure information collected while a closed-loop experiment are plotted in Figure 5. The hypoglycemia prediction alarm system advised rescue carbohydrates on a few occasions and prevented hypoglycemia during the 60-hour experiment in spite of daily exercise challenges. The adaptive controller worked well in keeping glucose concentration in the target range during most of the experiment. A critical contribution of this multivariable approach is observed on the last day during and after exercise (1 PM). Since the CGM values are high, a controller relying only on CGM values would have instructed the pump to infuse high levels of insulin. But since energy expenditure indicates physical activity, the controller takes into account the glucose-lowering effect of exercise and reduces insulin infusion rate to basal level to prevent a hypoglycemia episode.

Figure 4. Data collected in a clinical closed-loop control experiment. Glucose concentration from continuous glucose monitor (CGM), insulin infusion rate, and energy expenditure (EE). Hypoglycemia threshold: black dashed line. Vertical bars: black bar indicates regular meal, snack, or exercise (B, breakfast; L, lunch; D, dinner; S, snack; E, exercise), green (gray) bar indicates snack or juice based on hypoglycemia early alarm (J, orange juice) and magenta (gray) indicates sensor calibration (C).

Conclusions

Mathematical models are essential for many simulations, supervision, and control activities. Elucidation of the suitable type of model has a significant impact on the success of these activities. Two types of models, ABMs and time series models, illustrate the rich spectrum of modeling techniques and their contributions in developing powerful simulation and control algorithms that can aid in the advancement of the pharmaceutical industry and medical device development.

Acknowledgments

The authors are grateful to NSF for supporting of our ABM research (CBET-0731201, IIS-1125412) and to NIH (1DP3DK101075, 1DP3DK101077) and JDRF (17-2013-472) for supporting of our artificial pancreas research. The contributions of Dr. Arsun Artel and Dr. Hamid Mehdizadeh in the development of the earlier versions of the ABM of angiogenesis and the contributions of Dr. Meriyan Oruklu in the development of recursive time series modeling and adaptive control of glucose concentrations are acknowledged.

References

1. National Science Foundation. Simulation-Based Engineering Science: Revolutionizing Engineering Science through Simulation. NSF. 2006.
2. Casella G, Wu R, Wu AS, Weidman ST. Making Sense of Complexity: Summary of the Workshop on Dynamic Modeling of Complex Biomedical Systems. National Academies Press, Washington, DC. 2002.
3. Morreale P. Agents on the move. IEEE Spectr.1998;35(4):34-41.
4. Mehdizadeh H, Sumo S, et al. Three-dimensional modeling of angiogenesis in porous biomaterial scaffolds. Biomaterials. 2013;34(12):2875-2887.
5. Peirce SM. Computational and mathematical modeling of angiogenesis. Microcirculation. 2008;15(8):739-751.
6. Bayrak ES, Mehdizadeh H, et al. Agent-based modeling of osteogenic differentiation of mesenchymal stem cells in porous biomaterials. Engineering in Medicine and Biology Society (EMBC), 2014. 36th Annual International Conference of the IEEE, IEEE.
7. Zhang L, Chen LL, et al. “Multi-scale, multi-resolution brain cancer modeling.” Mathematics and Computers in Simulation. 2009;79(7): 2021-2035.
8. Thorne BC, Hayenga HN, et al. Toward a multi-scale computational model of arterial adaptation in hypertension: verification of a multi-cell agent based model. Front Physiol. 2011;2(20).
9. Kaul H, Cui Z, et al. A Multi-Paradigm Modeling Framework to Simulate Dynamic Reciprocity in a Bioreactor. PLoS ONE. 2013;8(3).
10. Artel A, Mehdizadeh H, Chiu YC, Brey EM, Cinar A. An agent-based model for the investigation of neovascularization within porous scaffolds. Tissue Eng Part A. 2011;17:2133-2141.
11. Eren-Oruklu M, Cinar A, Quinn L, Smith D. Adaptive control strategy for regulation of blood glucose levels in patients with type 1 diabetes. Journal of Process Control. 2009; 19(8):1333-1346.
12. Turksoy K, Quinn L, Littlejohn E., Cinar A. Multivariable adaptive identification and control for artificial pancreas systems. Biomedical Engineering, IEEE Transactions on. 2014;61(3):883-891.
13. Turksoy K, Bayrak ES, Quinn L, Littlejohn E, Rollins D, Cinar A. Hypoglycemia early alarm systems based on multivariable models. Industrial & engineering chemistry research. 2013;52(35):12329-12336.
14. Doyle FJ, Huyett LM, Lee JB, Zisser HC, Dassau E. Closed-loop artificial pancreas systems: engineering the algorithms. Diabetes care. 2014;37(5):1191-1197.
15. Bequette BW. Algorithms for a Closed-Loop Artificial Pancreas: The Case for Model Predictive Control. Journal of diabetes science and technology. 2013;7(6):1632-1643.
16. Turksoy K, Cinar A. Adaptive Control of Artificial Pancreas Systems-A Review. Journal of healthcare engineering. 2014;5(1):1-22.
17. Bakhtiani PA, Zhao LM, El Youssef J, Castle JR, Ward WK. A review of artificial pancreas technologies with an emphasis on bihormonal therapy. Diabetes, Obesity and Metabolism. 2013;15(12):1065-1070.
18. Cobelli C, Renard E, Kovatchev B. Artificial Pancreas: Past, Present, Future. Diabetes. 2011;60(11):2672-2682.
19. Turksoy K, Bayrak ES, Quinn L, Littlejohn E, Cinar A. Multivariable adaptive closed-loop control of an artificial pancreas without meal and activity announcement. Diabetes technology & therapeutics. 2013;15(5):386-400.
20. Turksoy K, Quinn LT, Littlejohn E, Cinar A. An Integrated Multivariable Artificial Pancreas Control System. Journal of diabetes science and technology. 2014;8(3):498-507.

Author Biographies

Ali Cinar is a professor of chemical engineering and biomedical engineering at Illinois Institute of Technology and the director of the Engineering Center for Diabetes Research and Education (ECDRE). His research specializations include diabetes and control of artificial pancreas systems; modeling of angiogenesis and tissue growth; modeling, supervision, and control with agent‐based systems; and multivariate statistics and artificial intelligence techniques. Dr. Cinar holds a PhD in chemical engineering from Texas A&M University. He is a Fellow of AIChE and member of ADA, BMES, IEEE, and TERMIS. He has published 2 books, over 180 technical papers in refereed journals and conference proceedings. The list of publications and detailed description of research interests is available at http://mypages.iit.edu/~cinar.

Elif Seyma Bayrak is a doctoral student in chemical engineering at Illinois Institute of Technology. Her research focuses on development of agent-based models for predicting interactive tissue growth and angiogenesis in three-dimensional scaffolds. She received a BS degree in chemical engineering at Hacettepe University, specializing in bioengineering.

Kamuran Turksoy is a doctoral student in biomedical engineering at Illinois Institute of Technology. His research focuses on system identification, adaptive control, development of multivariable control, predictive hypoglycemia alarm and fault detection and diagnosis systems for the artificial pancreas. He received a BS in electrical engineering at Yildiz Technical University, specializing in systems and control area.