5.1 Introduction

The need to control complex systems, both engineered and natural, pervades our lives. From the thermostat that controls the temperature in our homes to the software that controls flight characteristics of the space shuttle during landing, the vast majority of engineered systems have built-in control mechanisms. Being able to control certain biological systems is no less important. For instance, we control ecosystems for agriculture and wildlife management; we control different parts of the human body to treat and cure diseases such as hypertension, cancer, or heart disease. And we control microbes for the efficient production of a vast array of biomaterials. With control comes the requirement to carry it out in a fashion that is optimal with respect to a given objective. For instance, we want to devise a schedule for administering radiation therapy to cancer patients in a way that maximizes the number of cancer cells killed while minimizing side effects. We apply pesticides to fields in a way that minimizes environmental damage. And we aim to control the metabolism of engineered strains of microbes so they produce the maximal amount of biofuel. Thus, the need for optimal control is a problem we face everywhere. This chapter will focus on optimal control of biological systems.

The most common approach to optimal control is through the use of mathematical models, often consisting of one or more (ordinary or partial) differential equations. These equations model key features of the system to be controlled and include one or more variables that represent control inputs. The following example illustrates this approach; it is taken from [1], where more details can be found. The problem we want to focus on is the optimization of cancer chemotherapy, taking into account certain immunological activity. The two relevant variables are x, which represents the volume of the tumor to be treated, and y, which quantifies the density of so-called immunocompetent cells, capturing various types of T-cells activated during the immune reaction to the tumor. These two variables are governed by the system of ordinary differential equations,

Here, F is a function that represents the carrying capacity of the tumor, a concept adapted from ecology, in that it represents that ability of the environment to sustain the tumor. The Greek letters are constant parameters of the model. For instance, represents the rate of natural death of T-cells. For our purposes the exact model structure is not important, and details can be found in [1]. For instance, the term represents the observation that tumors suppress the activity of the immune system.

The control objective here is to reduce the tumor volume through the use of a chemotherapeutic agent. We will refer to this as the “treatment.” We implement this control with a term , where represents the control input, namely the amount of the chemotherapeutic agent administered. The factor x takes account of the assumption that the effect of chemotherapy on the tumor is proportional to the tumor volume. Thus, the first equation becomes

Here, represents the carrying capacity of the tumor, as above. The parameter captures the growth rate of the tumor, with dimension cells/day. The parameter denotes the rate at which cancer cells are eliminated as a result of activity by T-cells, so that the term in the equation models the beneficial effect of the immune reaction on the tumor volume. Finally, we assume that the elimination terms are proportional to the tumor volume. We therefore subtract a term in the x-equation.

What kind of drug regime is appropriate for a given cancer patient depends in part on the particular state of this patient’s disease. Our goal now is to optimize the influence of the control variable u, that is, the treatment, with respect to several factors. On the one hand we want to shrink the tumor through the administration of a maximal dose of the agent, and on the other hand we need to take into account the toxic side effects of the medication. Also, the treatment should be as short as possible. All this can be combined into the cost function

This equation represents the cost of applying treatment schedule u for a length of time T. The optimal control problem now is to find a control u that minimizes this cost function. Solving such control problems is an active area of research, and algorithms have been developed to find u[2]. See also [3] for optimal control problems in systems biology.

Modeling biological systems through differential equations has its limitations, however. In many cases, the processes involved might be fundamentally discrete rather than continuous. For instance, in the case of a predator–prey relationship between two species inside an ecosystem, both populations are comprised of discrete individuals that engage in typically binary discrete interactions. Thus, it is not immediately clear whether one can apply methods such as differential equations, which assume that the quantities modeled vary continuously. In molecular biology, when we study regulatory relationships between genes inside a cell, these relationships are based on the interactions of discrete molecules. Modeling such systems using differential equations is based on two assumptions: first, there are many individuals involved, so that we can view them collectively as a continuous quantity; second, that we are able to describe the individual interactions in a “global” manner, as a term in a set of differential equations, usually involving one or more global parameters. Sometimes, both of these assumptions are justified, such as for large populations of bacteria or large quantities of chemicals in a fermenter. But at other times, one or the other, or both, of these assumptions fail. For example, in a cell, there might only be two or three molecules of a particular protein present at any given time, so that its role in a regulatory network becomes discrete and stochastic, and cannot be accurately modeled with continuous functions. And, as another example, in an ecosystem with several species, continuous models sometimes cannot accurately capture extinction events when one species reaches a very low count.

For these reasons, and others, several other frameworks have been used to model biological systems. One of these consists of so-called “agent-based,” also called “individual-based” models. These types of model have a long history, going back to the 1940s and the work of John von Neumann [4].

Von Neumann was interested in the process of self-replication and aimed to construct a machine that could faithfully replicate itself. A theoretical instantiation of such a machine turned into the concept of a cellular automaton as a computational model. A very well-known example of a cellular automaton is John Conway’s Game of Life[5]. Since it includes many of the basic concepts of agent-based models, we describe it briefly.

The Game of Life takes place on a chess-board-like 2-dimensional grid. This grid can either be finite or extend infinitely in all directions, thereby yielding two different versions of the game. Each square, or cell, on the grid can be either black or white. Since the Game of Life is intended to simulate life, a cell is instead referred to as either alive (black) or dead (white). Each cell away from the boundary of the grid has eight neighbors on the grid that physically touch it, four with which it shares an edge, and four that touch only on the corners. Cells on the boundary have fewer neighbors. To make all cells uniform, one can make the assumption that a cell at an edge of the grid, but away from the corner, has as additional neighbors the corresponding cells on the opposite edge of the grid. Thus, we effectively “glue” opposite edges together, so that the grid is situated on a torus rather than in the plane. Now that all cells have eight neighbors, we establish a rule that determines the evolution of the cells on the grid by determining what the status of a given cell is, depending on the status of its eight neighbors. The rule is quite simple.

Thus, the rules are reminiscent of a population whose survival is affected by under- and overpopulation. If we now initialize this “Game” by assigning a “live” or “dead” status to each of the cells, then we can use this rule to evolve life on this grid by updating the status of all the cells in discrete time steps. The result is a vast array of different dynamics that can be observed, but is largely unpredictable from a particular initialization. There is a rich literature on this topic and many websites with Game of Life simulators. Before continuing, the reader is encouraged to try some of these and explore the question of predicting dynamic behavior from particular initializations.

General agent-based models have many features similar to this set-up. There is a collection of agents that are distributed across a spatial landscape. In the Game of Life, the spatial landscape is the grid, and there is one agent per cell. Each agent can be in one of a (typically finite) number of states, such as “dead” or “alive,” and has a set of rules attached to it that it uses to determine its state, based on the state of those other agents it interacts with at any given time. Beyond the Game of Life, in general agent-based models, agents are able to move around in the spatial landscape, and there are rules that determine the agents’ movement patterns. Typically, the rules are stochastic, rather than deterministic, and are governed by a collection of probabilities. For instance, agents might be predisposed to follow a certain gradient, representing, e.g., nutrient availability, but there is some chance that they might move in a different direction. While there are many other variations and features in particular agent-based models, these few basic features characterize the agent-based modeling framework.

These features are also sufficient to explain the basic differences between agent-based models and equation-based models, such as ordinary differential equations. Agent-based models lend themselves very well to a description of dynamical systems that arise from local interactions of many parts/agents, based only on local rules rather than on the configuration of the entire system at any given time. Also, it is very easy to represent a rich heterogeneous spatial environment that the agents navigate. Thus, the dynamics of the entire system, or its so-called global dynamics, “emerge” from these local interactions by applying the local rules repeatedly. In contrast, a system of differential equations, for instance, explicitly describes the global dynamics of the system up front. Furthermore, all the specifications for an agent-based model are intuitive, in the sense that they are direct computational representations of recognizable features in the actual system. This leads to models that are more faithful to the system to be modeled and that are more accessible to domain experts. With existing software for model building, they can even be built by domain experts directly, without the intervention of a modeler or mathematician. But, as always, there are no free lunches. These advantages come with some significant costs attached.

The reason that agent-based modeling became widespread only in the last 15 years or so is due to the fact that larger models with more features require extensive computational resources that were not available until recently. It is only now possible, barely, to build and analyze models that might have hundreds of millions of agents and tens of millions of spatial locations, with agents being very complex in terms of the states they can take on and the rules they follow. Even moderately complex models require high performance computing facilities for their analysis, which makes it difficult for individual researchers to use them. High computational complexity is compounded by the fact that there are very few mathematical tools available for the analysis of agent-based models. In essence, the only approach to model analysis is model simulation. That is, from given initializations of agent states and environmental parameters, one observes the time evolution of the system and hopes to observe patterns that might help one draw conclusions about, e.g., steady states of the model. Through choosing many initializations and doing many simulation runs from a given initialization in the case of a stochastic model, one aims to obtain global dynamic information about the model. Little else can be done because, in essence, the model consists of a computer program rather than a set of mathematical equations, and there are few things one can do with a computer program other than run it.

The lack of mathematical analysis tools extends in particular to the arena of optimal control. Existing approaches are heuristic in nature, based on domain knowledge. The goal of this chapter is to describe some of these approaches and to outline steps one can take to expand the repertoire of available tools to include techniques based on mathematics. The way to do this is to translate an agent-based model into a mathematical object, such as a system of equations of some sort, that makes it amenable to mathematical analysis tools. There are several possible ways in which to do this, and research is only in its early stages. Thus, the reader should see this chapter as a snapshot of an exciting research area that is evolving rapidly and providing rich opportunities for contributions at all levels. This chapter showcases one possible approach and the steps that have been accomplished on the road to developing mathematical tools for optimal control of agent-based models.

5.2 A First Example

Go to http://ccl.northwestern.edu/netlogo/models/RabbitsGrassWeeds. There you will find an agent-based model of rabbits in a field of grass and weeds. At each time step (or “tick”) the rabbits move in a random direction (they lose energy by moving). If there is grass at their location, they eat it and gain energy. If their energy level climbs above a certain threshold, they give birth (in this model, a new rabbit is spontaneously created at the location of the parent). Upon birth, the parents’ energy is halved, and the offspring is created with this halved energy level. Upon each tick, empty squares (or “patches”) have a certain fixed probability of grass re-generating. Weeds can also be introduced in order to further complicate the dynamics; their behavior is similar to that of the grass.

Near the top of the page you will find an option that allows you to run this model in your web browser. Spend some time reading through the description in order to get a feel for how this model works. Click setup and then go to run the simulation (press go again to stop the simulation). Note that you can speed up or slow down the model by using the slider at the top of the interface. Each time you wish to re-initialize the model and start over, you must click setup again.

In this model, each of the grid squares (henceforth referred to as patches) are agents, and each rabbit is an agent as well. Notice that the rabbits move around randomly, eating grass as they encounter it. Note too that the patches are colored green when grass is present, and black when it is absent. Notice the sliders on the left-hand side: these values can be changed as the simulation proceeds, or prior to starting a simulation. What do you expect to happen to the rabbit population if the grass growth is set to a very low rate? What if the grass grows quickly? What effect would altering the birth threshold have on the amount of grass present at each time step? Do you think it is possible for the rabbit population and grass level to stabilize over time? Play around with the sliders to determine if your predictions are correct. Note that you may need to let the simulation run for several hundred time steps (or “ticks” in Netlogo) in order to observe consistent dynamics. Do the populations stabilize? Do they oscillate?

5.3 Netlogo: An Introduction

Hopefully, the example in Section 5.2 has given you an idea of what an agent-based model is: a computer simulation wherein agents interact with their local environment (possibly including other agents) based on a set of rules. In this section we guide you through several exercises involving agent-based models using Netlogo [6], a software tool developed to work with agent-based models. Netlogo can be downloaded for free at http://ccl.northwestern.edu/netlogo/. There you will find detailed instructions on how to install it. Netlogo is its own programming language, so named because it is a variation of the Logo language. While many key features are unique to Netlogo, users familiar with Logo will likely note similarities. Note too that Netlogo is continually updated and newer versions released. As of the time of this writing, the latest version is Netlogo 5.0. All files and exercises associated with this chapter will be conducted using this version; some adjustment may be necessary for newer (or older) versions. While there are many other software platforms available and in use for agent-based modeling, for the remainder of this chapter we will use Netlogo to introduce and examine agent-based models. It has the convenient advantage of providing the user with an intuitive graphical interface, which we will use to aid our understanding of the models we will examine. In addition, the standard installation comes with a library of models, all of which are open source. Thus we may build models from scratch or we may choose to alter an existing model. You may wish to go through the tutorials found at http://ccl.northwestern.edu/netlogo/docs/. A complete dictionary of programming terms can be found at http://ccl.northwestern.edu/netlogo/docs/dictionary.html. All exercises in this section refer to the “Rabbits Grass Weeds” model introduced in Section 5.2. It should be noted that for Exercises 5.1 and 5.2, the web-based version of the model can be used (via http://ccl.northwestern.edu/netlogo/models/RabbitsGrassWeeds), but Exercises 5.3 and 5.4 require installation of the Netlogo software.

5.4 An Introduction to Agent-Based Models

Mathematical modeling is a method of encoding relationships and interactions in a natural or engineered system into a formalized system; the models can then be studied and analyzed using a variety of mathematical approaches. Models allow researchers from a wide variety of disciplines to examine systems and their emergent behavior. For example, a model may be used in order to make predictions about the future behavior of a system, or it may be used to solve a complicated problem explicitly. In developing the best drug to administer to treat a particular disease, researchers in the past have relied on trial-and-error methods via repeated experimentation on live subjects. With an accurate model of the subject, however, such experiments can often be faithfully reproduced in silico; that is, they can be run on a computer and analyzed with far less preparation and expense. This highlights another benefit to mathematical modeling: given an accurate model of the system in question, a year’s worth of real time can be simulated in several minutes.

Experimentation remains part of the process - a model will always have a limited ability to predict and must be correlated against empirical data in order to ensure that the models are indeed faithful simulations of the actual physical system. Thus, models are best used to inform in vivo experimentation, which in turn produces results that can be used to calibrate the model further.

Agent-based models are a class of computer models in which entities (referred to as agents) interact with each other and or their local environment. Formally:

Systems (such as the human immune system) are increasingly being implemented in the form of agent-based models (with individual cells as the agents, of which there may be many types) as more and more research involves the use of in silico simulation to study the properties of this and other similarly complex systems.

Agent-based models have the advantage of being well suited for modeling many different types of systems. They have been used to study social interactions among individuals, the spread of disease through populations, scheduling and efficiency of factory processes, how cells react to drug treatments, and many other systems. It is perhaps worth noting that many of the current issues in the scientific community are interdisciplinary. Finding a cure for cancer will involve geneticists, biologists, mathematicians, chemists, and perhaps many other specialists. Agent-based models have an intuitive formulation and can often be examined via a graphical interface. Thus they are a natural tool for promoting interdisciplinary research, as the mathematics underlying the models is hidden in the programming; in other words, it is possible for biologists, chemists, and other researchers to make use of agent-based models without a full background in the mathematics that are involved in creation and analysis of the model. At the same time, the mathematical structure remains and can be explored concurrently by mathematicians.

Even within the scope of mathematical modeling, there are several noteworthy advantages to agent-based modeling. One such advantage is that they are effective for modeling systems wherein many agents follow the same set of rules (e.g., rabbit populations or blood cell types). In particular, models containing spatial heterogeneity can be effectively represented via agent-based models while perhaps not so easily using, for example, ordinary or partial differential equations (ODEs). Altering the spatial dynamics of an agent-based model consists of changing several lines of code or less, while such changes can be difficult (and perhaps even impossible) to implement in an ODE model.

While agent-based models provide a convenient and natural setting for studying complex systems, there are several issues within the research community that are currently unresolved. A 2006 paper [8] written by a large group of researchers identified two main obstacles: the first is the lack of standardization of the description of agent-based models, and the second is a lack of rigorous mathematical formulation of the system itself. Descriptions of agent-based models can vary in different settings, and as of this writing there is no standard definition that is universally agreed upon. Some models are developed to simulate physical processes, and others are developed in the framework of graph theory. In some cases models are developed in order to study only a certain aspect of the given system; thus the model may have more variables pertaining to certain processes than others. A research article may begin with the statement of an objective, or it may begin with a description of the model itself. In some cases, the rules that govern the updating of the agents’ state variables are deterministic, and in other cases they are stochastic. All of these issues would be resolved if there was a standard protocol for describing agent-based models (indeed, one such protocol is proposed in [8]). The need for an agreed-upon structure to be followed is perhaps most clearly felt in the model’s presentation. In particular, the layout of the model presentation ought to be standardized so that a reader immediately knows where to look in the description in order to learn what the model describes and what its rules are. In the current literature models are presented in myriad forms, and the descriptions of the agents, environments, and rules come in no particular order. Thus a reader is required to scour the paper for pertinent details that might otherwise be presented in some standard way.

The second major issue concerns the lack of rigor in the formulation of the model. In many cases, the description of a model is given in several paragraphs, describing in some imprecise manner what the agents are and how they interact with each other. In fact, in order for agent-based models to be implemented on computer software, there are intricate rules embedded in the programming of the models. These rules and equations are often glossed over if not entirely omitted; if they are given, it is typically through a verbal description rather than a strict mathematical formulation. Looking back on the example in Section 5.2, the precise way in which the rabbits move is not described. For example, it is unclear from a simple description whether the rabbits can move diagonally, or whether the distance they move at each time step is variable or constant. In fact, even running the model does not immediately answer this question—it is only made explicitly clear by thorough examination of the computer code. In order to describe such a model rigorously, it is necessary that this information (and other similarly imprecise descriptions) be presented clearly and unambiguously.

In addition to these issues, an author may spend several paragraphs explaining how the model was formulated that could just have easily been given in one or two equations. Short and precise definitions can save time for the reader and also make the author’s intentions explicit. One proposition to overcome this lack of formalization is presented in [9]; this work proposes a rigorous mathematical representation for agent-based models as polynomial dynamical systems. A further description of such systems can be found in [10] and is described in Section 5.8.

5.5 Optimization and Optimal Control

Agent-based models are typically created to simulate some real-world process in order to aid investigation. Once a model has been created, a natural next step is to ask: what are we to do with it? In this section we give a brief overview of optimal control and optimization, and introduce these ideas as they apply to agent-based models.

Optimization is the process of finding the best solution with respect to a particular goal. For example, suppose we have a model of the immune system battling a bacterial infection, and we wish to study the effects of certain drugs on this battle. We may wish to find out which drug does the least amount of tissue damage while curing the infection—this is an optimization problem, because we are searching for the best drug with respect to the stated goal. On the other hand, perhaps our goal is to cure the infection in the shortest time possible, regardless of the tissue damage caused. This would be another optimization problem. It is likely, then, that the solution to the optimization problem depends on the optimization goal. In this scenario, the drug which cures the infection most quickly may consequently do more tissue damage than other drugs (though not necessarily); on the other hand, the drug that causes the least tissue damage might not cure the infection in the shortest possible time. In this example, optimization is a process of minimization: in one case we wish to minimize tissue damage and in the other we wish to minimize the healing time. However, it is also possible that the solution of an optimization problem is a process of maximization: given a model of an immune system fighting a fatal illness, what treatment will enable the patient to survive the longest? In fact, it may be that the goal of the optimization problem is to minimize one value while maximizing another. For example, our optimization goal may be to minimize tissue damage and maximize the expected lifespan of the patient, given that we can only administer a particular drug a fixed number of times. In general, however, optimization is the process of finding the best solution, depending on the objective.

Typically, optimization is a complicated process, and becomes even more so when realistic constraints are put in place. Through the use of agent-based models it may be possible to obtain a solution to an optimization problem that is not feasible in actuality: for example, the solution may exceed monetary limitations, may require actions that are not permitted by health care regulations, or may require interaction with the patient in an impractical way (e.g., if the solution calls for treatment every hour for 100 consecutive hours, or for of the population to receive vaccination). Nevertheless, optimization remains a natural means of applying mathematics to solve real-world problems. Such problems are often framed as questions of optimization: what is the best outcome that can be produced based on the properties of the model?

Optimal controlis a slightly different notion. In an optimal control problem, we ask the following question: given a particular state that we hope to reach, what is the most efficient way of reaching that state? In other words, we know the state of the system that we hope to reach, and the control problem is to find the solution that steers the system to that state in the best manner. Using a model of the immune system fighting an infection as an example, we may formulate an optimal control problem as follows: given that we wish to eliminate the number of infected cells in one month, what is the best drug treatment schedule we can devise? To summarize: in an optimal control problem we know the state we hope to reach and we search for the best way to reach it, whereas in an optimization problem, a goal is stated and we compare solutions to see which solution maximizes or minimizes the stated goal of the problem.

We now present an example explaining several terms related to optimization and optimal control that will be used in Sections 5.6 and 5.7. We then conclude this section with definitions of said terms. We do this in order to standardize and clarify the terminology within this chapter. Note that while many of these terms are also used in optimal control for continuous systems, there might be subtle differences in their meanings when applied to agent-based models. Furthermore, as this chapter is meant to provide an introduction to the topic, more formal definitions are outside the scope of this text. For a more formal treatment see [11].

In this case, there will be many variables: the number of healthy cells, the number of cancer cells, the rate at which cancer cells grow, the rate at which healthy cells regenerate, the expected lifespan of the patient, the frequency with which the drug is administered, the type of drug, and so on. Some of these variables will have fixed values: for example, the rate at which healthy cells regenerate (during intervals when no treatment is administered) can be determined through experimental measurements, and in fact this rate helps to define the model itself. Such variables are referred to as model parameters—they are a part of the specification of the model. The repeated interactions of the entities in the model, such as immune cells or rabbits, are referred to as the model dynamics. Note that we will only have direct control over some of the variables, and others we will simply measure by observation. For example, during each day of the simulated treatment we can decide whether or not to administer the drug; thus we have direct control over the value of this variable. We refer to these as control variables, because we have direct control over their values and they exercise control over certain aspects of the model. On the other hand, we might not be able to control over other aspects, such as the number of white blood cells present at the site of an infection. We refer to variables of this type as state variables because they help to describe the state of the system at any given time.

In this example, our only control variable is a binary decision whether or not to administer the drug on a given day, thus there are only two possible values for this control variable at each time step. For any given day, we may represent that the drug was administered that day by assigning that variable the value 1, and represent that the drug was not administered by assigning that variable the value 0. In this situation, then, a treatment schedule is a vector of length 365, of which each entry is either 0 or 1. Thus there are a total of possible treatment schedules. However, it is likely that not all possible treatment schedules will reduce the number of infected cells below the fixed level. The treatment schedules that do achieve this are solutions to the optimal control problem. The solution space is the set of all such solutions. Note that this space is entirely separate from the state space, which consists of all possible states a patient may be in throughout any simulation of the model.

Recall that we are trying to find the best solution among many candidates. In order to do this, we must have some way of ranking individual solutions. We do this by introducing a cost function. We are attempting to steer the system to a particular state, and each solution achieves that goal at a certain cost. We wish to determine the best solution (i.e., the optimal control); thus we wish to find the solution that minimizes the cost function. Earlier we noted that in the immune system example the best solution was the solution that involved the least number of treatment days throughout the treatment schedule; thus, in this example a reasonable cost function might involve the sum of the entries in the treatment schedule (i.e., the number of days on which the drug was administered). For the purposes of this chapter, the goal of the optimization and the optimal control problems will always be to minimize the associated cost function.

In our example, we wish to minimize T, the total number of days the drug is taken, and maximize I, the number of healthy cells. Then a reasonable cost function for a treatment schedule S which contains days on which the drug is administered and results in an expected value of healthy cells might be . Here we want to maximize over all treatment schedules S; this is the same as minimizing (in this way we can always formulate our problem so that the goal is to minimize the cost function). However, it is perhaps more realistic that one of our optimization goals has a higher priority than the other. For example, we may suppose that maximizing the number of healthy cells is more important than minimizing the number of treatment days. In that case, we introduce weighting coefficients to alter the cost function. For example, it may become . This has the effect of increasing the relative importance of the healthy cell count when evaluating the associated cost of a particular treatment schedule.

In order to use our model to obtain results, we simulate the model a number of times and make observations on the results. We might, for example, administer a treatment every day. This particular solution is represented as a string of 365 consecutive 1’s; we implement this in the agent-based model and count how many healthy cells are there after one year of simulated time. We then evaluate the cost function based on this information. Given that agent-based models are often stochastic in nature (meaning that while the rules are fixed, the model dynamics depend on many random variables), results obtained from one simulation are typically not reliable. We eliminate such variation by simulating a given control schedule many times and using averages for evaluation in the cost function.

Of all the possible solutions, there must be one that is better (or at least no worse) than every other solution. Strictly speaking, such a solution is referred to as an optimal solution. However, in many cases we cannot simulate all possible solutions because of the monumental amount of computing time involved. Thus, for the purposes of this chapter, we refer to the best solution we are able to find as the optimal solution, with the caveat that there may indeed be a better solution elsewhere in the solution space.

We now define the following terms, each of which have been illustrated in this example.

Use the modified rabbit and grass model file RabbitsGrass.nlogo at http://admg.vbi.vt.edu/software/rabbitsgrass-netlogo-files/to answer the questions in this section. Once you have opened the file, you will find a description of the interface settings under the Info tab.

5.6 Scaling and Aggregation

Our investigation of agent-based models and how to formulate control problems for them motivates this section and the next. Searching for an optimal solution in such models often requires running many thousands of simulations, thus performing such simulations as quickly as possible is a primary concern. In this section, we discuss means of reducing the run time and complexity of agent-based models via scaling and aggregation.

Scaling is a method of shrinking the size of an agent-based model in order to improve run time. In Exercise 5.1, we explored the “Rabbits Grass Weeds” model. The dimensions of this model are patches for a total of 1935 patches. The default number of rabbits is 150. In an attempt to scale this model, we may reduce the dimensions (and correspondingly, the initial number of rabbits). If we change the dimensions to for a total of 625 patches, we may choose the initial number of rabbits to be . Our hope, then, is that the pertinent model dynamics remain the same at this reduced size, since in that case we can run all subsequent trials at this reduced size for a substantial decrease in run time.

The first question we need to answer is exactly what it means for the model dynamics to remain the same, and how we can verify that this is indeed the case? Since we are using the model with a specific control objective in mind, and the value of a solution relies on the associated cost function, this cost function will help us determine how to quantify whether our scaled model retains the pertinent dynamics. In Section 5.5, we used the “Rabbits and Grass” model to determine control schedules for poisoning the rabbits. In that case, our cost function relied on two parameters: the rabbit population, and the number of days in which poison was used. In particular, the day-by-day grass levels are not relevant to the cost function, nor are the energy levels of the rabbits; thus these parameters need not be preserved in the scaled model. In the scaled model, the number of days in which poison is used is unchanged, so this variable is also unaffected by scaling. Then the only parameter we need to preserve is the rabbit population. To be more specific, we need to answer the following question: what are the smallest dimensions we can use in the model (by scaling down the number of rabbits accordingly) so that the average day-by-day rabbit population dynamics in the scaled model follow the same pattern as those dynamics in the original model? We answer this question by simulating many runs at each size and keeping track of the population levels for each. We can then calculate the correlation using any of a number of statistical methods; one such method is now presented.

5.7 A Heuristic Approach

In previous sections we have discussed optimal control as it relates to agent-based models, and techniques for improving the run time of such models. In theory, the process of finding the optimal solution requires only that we evaluate the entire solution space and choose the solution that minimizes the cost function. In practice, however, this remains unfeasible. For one thing, it is quite possible that the solution space is infinite; and even if it is finite, it is often the case that there are far too many solutions to possibly implement them all (in terms of realistic run times), even in models that have been subjected to scaling and aggregation techniques.

For example, a discrete model of the immune system might involve many millions of cells. Each will have its own rules describing its interaction with other cells, and may also have many variables attributed to it. Analytic methods fail in these cases, and exhaustive enumeration by computer is unfeasible due to the limitations of computer processing speeds. In fact, as the processing power of computers grows, so too does the possibility for creating models that increase in complexity, so that it is impossible to be certain that computers will ever be able to “catch up” to the complexity of the models.

For this reason, most of the optimization and optimal control methods currently used in discrete modeling (and with agent-based models in particular) are in the form of heuristic algorithms. An algorithm is a formal set of rules used in order to obtain a solution to a problem. A heuristic algorithm is a specific type of algorithm that conducts a “smart” search of the solution space. It may make use of certain aspects of the problem to avoid having to search through every possible solution, or it may update its search method based on the input it receives from previously found solutions. These algorithms begin with a choice of control input values (i.e., full or partial solutions) and use various methods to refine the values so as to decrease the value of the associated cost function until no better solution can be found. While these methods do not guarantee an optimal solution, they do provide an opportunity for analysis. In particular, if constructed and executed correctly, they are a vast improvement upon random searches of the solution space. One advantage to the heuristic approach is that it can be used with virtually any model. Since heuristic algorithms rely on results from simulation, it is not necessary to transform the model in any way in order to implement them. They provide a “brute force” technique that, while not always optimal, can always be used when other methods may fail. In particular, heuristic algorithms provide a baseline for results obtained via other, more rigorous methods. While a mathematician may be concerned with developing an algorithmic theory for explicitly finding a true optimal solution, in practice scientists and other researchers are often satisfied with a solution that is better than any previously known. As such, the heuristic approach to optimal control is a valuable option.

Of course, every heuristic algorithm has associated risks and advantages, and the particular algorithm that is best suited for a given model or problem will depend on these. As an introduction to the heuristic approach, we describe a so-called genetic algorithm in some detail, in the context of the “Rabbits and Grass” model. We conclude this section with a list of other heuristic algorithms.

5.8 Mathematical Framework for Representing Agent-Based Models

As addressed earlier, many agent-based models are too complex to find global optimal control purely by simulation. One option one might pursue in this case is to translate agent-based models into a different type of mathematical object for which more mathematical tools are available. This is similar to the approach pioneered by Descartes and his introduction of a coordinate system. In the plane, for instance, a Cartesian coordinate system allows us to represent lines by linear equations, circles by quadratic equations, etc. If we now want to find out whether two lines in the plane intersect we can either carry out an inspection of the two lines as geometric objects or, alternatively, we can determine whether the resulting system of linear equations has a solution. Translating agent-based models into a formal mathematical framework then opens up the possibilities to use theory from mathematics to analyze them from a different angle rather than just with simulations. One framework that has been studied is that of so-called polynomial dynamical systems (PDSs).

As in our running examples of rabbits and grass, agents in the models we consider here take on a finite number of states, such as a rabbit’s position on the grid or the amount of grass on a patch. This finite set of states can be viewed as the analog of the points that make up a line in the plane, to continue with the analogy of the Cartesian coordinate system. Applying the rules of the model in order to update the states of all the agents can be thought of as a function that maps the collection of discrete states to itself. Since the state sets do not carry any kind of structure we are simply dealing with a function from a finite set to itself, about which we can say little. However, let us now introduce the analog of a Cartesian coordinate system into this set in a way that provides extra tools for us, analogous to the ability to use algebra to decide whether two lines in the plane intersect. We do this by way of an example.

Suppose that our set of states equals , which we can represent as {0, 1, 2}. We can easily turn this set into a number system by using so-called arithmetic modulo 3. That is, we can add and multiply the elements of this set according to rules that are identical to those used to multiply real numbers. Notice that the numbers in the set are all the remainders one can obtain under division by 3. Now we add and multiply “modulo 3,” that is, 2 + 2 = 1 and as well, that is the remainder of 2 + 2 = 4 under division by 3 is 1. The number system we obtain in this way is an example of a “finite field,” which we will use for the remainder of the chapter. We denote the finite field with three elements by , or more generally, the finite field with p elements, where p is prime, by .

Given rules from an ABM that describe how the agents or variables change their states, we can find for every agent a polynomial that describes the given rule using elements of the finite field such as {0, 1, 2} rather than {low, medium, high}. By constructing a polynomial for every single agent, we obtain a set of polynomials that describes the same behavior as the original ABM. The set is called a polynomial dynamical system.

Here, denotes the Cartesian product with n copies of k. denotes the ring of polynomials in variables and coefficients in k. For the function F is evaluated by

The PDS F maps points in to other points in , and iteration of F results in a dynamical system. F has one fixed point, or steady state, namely (1, 0, 2) because , and note that all calculation are done “modulo 3” because F is over . Furthermore, F has a limit cycle or oscillation of length five: applying F repeatedly to (0, 2, 1) yields

With software packages such as ADAM [25], one can compute the dynamics of a PDS and visualize it.

Most agent-based models can be translated into a polynomial dynamical system of this type. Each variable represents an agent, a patch, or some other entity in the model. Each element in the finite field k represents a different state or condition of the variables. The corresponding polynomials encode the behavior of the agent or patch as described in the model. In the rabbit and grass model, one could assign a variable for each patch. The values of the finite field then represent how many rabbits and how much grass is currently on the patch. The polynomials determine the number of rabbits and amount of grass on patch i at the next iteration, given the current number of rabbits and grass. We will explain in more detail how to derive the polynomials for an agent-based model in Section 5.9, however their existence is assured by the following result.

Any such mapping over a finite field can be described by a unique polynomial. Using Lagrange interpolation, we can easily determine the polynomial. Let be any function on k. Then

(5.1)

is the unique polynomial that defines the same mapping as f.

Then the polynomial g that defines the same mapping as f is constructed as follows:

Converting an agent-based model into a polynomial dynamical system provides us with a conceptual advantage, since rather than being limited to working with a computer simulation as our only means of analysis, methods and theory from abstract algebra and algebraic geometry can be used. For example, one might be interested in the steady states of a model, i.e., all the configurations of the system that do not change over time. Written as a polynomial dynamical system , these states are exactly the solutions to , the n-dimensional system of equations in [25]. Equivalently, the solutions are the points in the variety of the ideal generated by the polynomials . For an introduction to varieties, we recommend [27].

When the polynomials describing the biological system have a special structure, other analysis methods are available. For example, for conjunctive networks, i.e., a PDS over where each polynomials is a monomial, the dynamics can be completely inferred by looking at the dependency graph or wiring diagram of the PDS [28]. We can determine the fixed point and limit cycle structure of this PDS without using any simulation.

We have mentioned Conway’s Game of Life, a cellular automaton as a special case of agent-based models earlier in this chapter. Cells are agents that die or come to life based on a rule including their eight neighbors [5]. The Game of Life can be translated into a polynomial dynamical system. Each variable represents a cell on the grid. Each polynomial depends on ’s eight neighbors and itself and describes whether will be dead (0) or alive (1) at the next iteration given the values of its neighboring cells. The fixed points of this system correspond to still lives such as blocks and beehives, two cycles to oscillators, e.g., blinkers and beacons. Working with the mathematical representation of the game, one can for example study still lives by using concepts from invariant theory using the symmetry of the rules as group actions.

In the next section, we provide some polynomials for common rules in agent-based models. By providing such rules, we hope to simplify the process of translating an ABM to a PDS. The polynomials can be used as given or as a starting point to construct functions that represent more complex behavior.

5.9 Translating Agent-Based Models into Polynomial Dynamical Systems

This section is authored jointly by Franziska Hinkelmann, Matt Oremland, Hussein Al-Asadi, Atsya Kumano, Laurel Ohm, Alice Toms, Reinhard Laubenbacher

Discrete models, including agent-based models, are important tools for modeling biological systems, but model complexity may hinder complete analysis. Representation as a PDS provides a framework for efficient analysis using theory from abstract algebra. In this section, we provide polynomials that describe common agent interactions. In the previous section we described how to interpolate agent behavior and to generate the appropriate polynomial. However, for a variable with many different states, this method of interpolation results in long and complex polynomials that are difficult to expand, simplify, or alter. Thus we provide some general “shortcut rules” for constructing polynomials that describe key agent and patch interactions present in many ABMs. Each of the following polynomials exists in the finite field .

Since we are particularly interested in ABMs describing complex biological systems, we use the term concentration to describe the states of a patch variable (for example, concentration of white blood cells on a patch). In this chapter we describe several polynomials that describe both basic movement, and movement according to the state of the neighboring patches.

5.10 Summary

Agent-based models provide a very intuitive and convenient way to model a variety of phenomena in biology. The price we pay for these features is that the models are not explicitly mathematical, so that we lack mathematical tools for model analysis. For instance, many of these phenomena are connected to optimization and optimal control problems, as pointed out in this chapter, but no systematic methods are available for agent-based models to solve these. We have attempted here to do two things. Firstly, we described so-called heuristic local search methods, such as genetic algorithms, which can be applied directly to agent-based models. And we described a way in which one can translate an agent-based model into a mathematical object, in this case a polynomial dynamical system over a finite field. Many computational and theoretical tools are available for such systems. For instance, to compute the steady states of a polynomial system , one can solve the system of polynomial equations

There are several computer algebra systems available to solve such problems. To compute the steady states of an agent-based model, on the other hand, one is limited to extensive simulation, without any guarantee of having found all possible steady states.

The chapter provides a snapshot of ongoing research in the field. The approach via polynomial dynamical systems, for instance, is very promising, but still lacks appropriate algorithms that scale to larger models. In addition to searching for such algorithms, further research in model reduction is taking place, as outlined earlier in the chapter. At the same time, other mathematical frameworks, such as ordinary or partial differential equations and Markov models are being explored for this purpose. Much work remains to be done but, in the end, a combination of better algorithms, improvements in hardware, and dimension reduction methods is likely to provide for us a tool kit that will allow the solution of realistic large-scale optimization and optimal control problems in ecology, biomedicine, and other fields related to the life sciences.

5.11 Supplementary Materials

All supplementary files and/or computer code associated with this article can be found from the volume’s website http://booksite.elsevier.com/9780124157804

References

1. Ledzewicz U, Naghnaeian M, Schattler H. Optimal response to chemotherapy for a mathematical model of tumor-immune dynamics. J Math Biol. 2012;64:557–577.

2. Sontag ED. Mathematical control theory. 2nd ed. New York, NY: Springer; 1998.

3. Iglesias PA, Ingalls BP, eds. Control theory and systems biology. Cambridge, MA: The MIT Press; 2009.

4. Von Neumann J. The general and logical theory of automata. In: New York: John Wiley & Sons; 1951;Jeffres LA, ed. Cerebral mechanisms of behavior – the Hixon symposium. vols. 1–31.

5. Gardner M. The fantastic combinations of John Conway’s new solitaire game life. Sci Am. 1970;223:120–123.

6. Wilensky U. NetLogo Center for connected learning and computer-based modeling. Evanston, IL: Northwestern University; 2009; In: http://ccl.northwestern.edu/netlogo/; 2009.

7. Laubenbacher R, Jarrah AS, Mortveit H, Ravi SS. Mathematical formalism for agent-based modeling. In: Springer 2009;Meyers RA, ed. Encyclopedia of complexity and systems science. vols. 160–176.

8. Grimm Volker, Berger Uta, Bastiansen Finn, et al. A standard protocol for describing individual-based and agent-based models. Ecol Model. 2006;198:115–126.

9. Hinkelmann F, Murrugarra D, Jarrah AS, Laubenbacher R. A mathematical framework for agent based models of complex biological networks. Bull Math Biol 2010.

10. Alan Veliz-Cuba, Jarrah Abdul Salam, Laubenbacher Reinhard. Polynomial algebra of discrete models in systems biology. Bioinformatics. 2010;26:1637–1643.

11. Kirk Donald E. Optimal control theory: an introduction. Englewood Cliffs, New Jersey: Prentice-Hall Inc; 1970.

12. Goldberg DE. Genetic algorithms in search, optimization, and machine learning Reading. MA: Addison-Wesley; 1989.

13. Holland John H. Adaptation in natural and artificial systems. Ann Arbor, Mich.: University of Michigan Press; 1975; [An introductory analysis with applications to biology, control, and, artificial intelligence 51].

14. Koza John R. Evolution and co-evolution of computer programs to control independently-acting agents. In: MIT Press 1991; Proceedings of the First International Conference on Simulation of Adaptive Behavior. vols. 366–375.

15. Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by simulated annealing. Science. 1983;220:671–680.

16. Pennisi M, Catanuto R, Pappalardo F, Motta S. Optimal vaccination schedules using simulated annealing Bioinformatics. 2008;24:1740–1742.

17. Fred Glover. Tabu search–part I. Informs J Comput. 1989;1:190–206.

18. Fred Glover. Tabu search–part II. Informs J Comput. 1990;2:4–32.

19. Glover Fred. Tabu search: a tutorial. Interfaces. 1990;20:74–94.

20. Dorigo M, Maniezzo V, Colorni A. Ant system: optimization by a colony of cooperating agents Systems. IEEE Trans Syst Man Cybern B Syst Humans. 1996;26:29–41.

21. Dorigo Marco, Caro Gianni Di, Gambardella Luca M. Ant algorithms for discrete optimization. Artif Life. 1999;5:137–172.

22. Joslin David E, Clements David P. Squeaky wheel optimization. J Artif Intell Res. 1999;10:353–373 (electronic).

23. Li Jingpeng, Parkes Andrew J, Burke Edmund K. Evolutionary squeaky wheel optimization: a new analysis framework evolutionary computation; 2011. published online.

24. Stigler B, Jarrah A, Stillman M, Laubenbacher R. Reverse engineering of dynamic networks. Ann NY Acad Sci. 2007;1115:168–177.

25. Hinkelmann Franziska, Brandon Madison, Guang Bonny, et al. ADAM: analysis of discrete models of biological systems using computer algebra BMC. Bioinformatics. 2011;12:295.

26. Laubenbacher Reinhard, Stigler Brandilyn. A computational algebra approach to the reverse engineering of gene regulatory networks. J Theor Biol. 2004;229:523–537.

27. Cox David A, Little John, O’Shea Donal. Ideals, varieties, and algorithms: an introduction to computational algebraic geometry and commutative algebra, 3/e (Undergraduate Texts in Mathematics). Secaucus, NJ, USA: Springer-Verlag New York Inc; 2007.

28. Jarrah Abdul, Laubenbacher Reinhard, Veliz-Cuba Alan. The dynamics of conjunctive and disjunctive boolean network models. Bull Math Biol. 2010;72:1425–1447. doi 10.1007/s11538-010- 9501-z.