DATA, MACHINE LEARNING, AND POLICY INFORMED AGENT-BASED
MODELING

By

David J. Butts

A DISSERTATION

Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

Computational Mathematics, Science, and Engineering – Doctor of Philosophy

2024

ABSTRACT

Agent-based models (ABMs) examine emergent phenomena that arise from individual
agent rules. This work extends the basic ABM paradigm in three key areas: data integra-
tion, evaluation of policies, and incorporating machine learning techniques. The dissertation
investigates how data-driven approaches can enhance the accuracy of ABMs, explores the
practical applications of ABMs in developing policies for real-world issues, and examines the
fusion of machine learning with ABMs to optimize model design and functionality.

The dissertation begins by establishing the background and fundamental principles of
agent-based models, highlighting their evolution from simple cellular automata to intricate
systems encapsulating decision-making and adaptive behavior. It examines the role of these
models in simulating dynamic interactions within systems, especially in scenarios where tra-
ditional methods may fall short of capturing the complexities of agent interactions. Following
the introduction of key concepts, a series of projects that function in pairs demonstrate the
versatility of ABMs and address the three key areas discussed above.

The ﬁrst pair addresses data integration of GPS deer movement data into a generalized
Langevin model and its use in uncertainty quantiﬁcation of disease spread. Exploratory
data analysis revealed a discernible non-parametric trend in the GPS data with non-Gaussian
statistics. This analysis led to a model that is consistent with the observed data. Subsequent
incorporation of chronic wasting disease (CWD) and population dynamics were used to
forecast the prevalence of CWD. This extended model was analyzed with a global sensitivity
analysis that tied variance in disease prevalence to variance in the parameters of the model,
providing predictions of future prevalence of the disease.

The second pair examines policy evaluation, speciﬁcally strategies for mitigating disin-
formation in social networks. Multiple strategies were evaluated on various topologically
diverse networks that led to policy recommendations. Simulations on these graphs revealed
challenges associated with large network simulations, particularly in computational cost and
inﬂuence of network topologies. These challenges led to a method to miniaturize real social
networks while preserving key attributes, enabling more eﬃcient and realistic simulations to
run on artiﬁcial social networks.

The ﬁnal pair investigates the possibility of inverting the ABM paradigm to instead have
agents learn their own rules through environmental interactions. Reinforcement learning
was applied to a model of conﬂict based on capture the ﬂag, where an agent learned in
progressively diﬃcult competitions. The emergence of deterrence was explored through
adding asymmetries between competing teams, and diﬀerential equation-based models were
created to help interpret results.

To Guido van Rossum, Leslie Lamport, and Linus Torvalds.
Without you, creating this dissertation would have literally been impossible.

iii

ACKNOWLEDGEMENTS

Throughout my time spent in graduate school, I have received a large amount of support
from many people. I am sure I have not listed everyone who has helped me along the way,
but that does not mean I do not appreciate your support.

First and foremost, I want to thank to my advisor Michael Murillo. Your guidance,
I
mentorship, and friendship has shaped me into a better scientist and a better person.
will greatly miss the many conversations had over wine and old fashions, and time spend
kayaking and walking around campus.
I also want to thank the remaining members of
my PhD committee. Arika Ligmann-Zielinska, I really enjoyed your class in agent-based
modeling, and your continued assistance with the sensitivity analyses in my projects. John
Luginsland, it is always fun to have conversation with you, and I hope to be as humorous
as you when I have a “real job”. More seriously, thank you for your support when I was
searching for jobs. Yuying Xie, you taught what I can only describe as the hardest class I
have ever taken. However, I appreciate it because not only did I learn a lot during your class,
but I gained conﬁdence in my abilities to solve diﬃcult problems. To each of you on my
committee, thank you again for your guidance and support throughout my time in graduate
school.

I also want to thank the Murillo Group’s current and past members. Speciﬁcally, Zach
Johnson, Jorge Martinez-Ortiz, Luciano Silvestri, Jannik Eisenlohr, Chris Gerlach, Thomas
Chuna, Luke Stanek, and close collaborators Liam Stanton and Jeﬀ Haack. You have all
helped me improve my research and presentation skills through countless meetings and prac-
tices. I also really enjoyed our many happy hours, including the ones that occurred virtually
throughout the pandemic. It was great being part of a research group that eﬀortlessly dou-
bled as a group of friends. Outside of the Murillo group, I want to thank Cole Stewart who
volunteered to read and provide feedback for many of the chapters in this dissertation.

I would like to thank the the many collaborators here at Michigan State and at Los
Alamos National Lab who have been a part of the research presented in this dissertation
and research I have done outside what is discussed here. I would deﬁnitely not be where I
am now without you, and it has been great learning about your areas of expertise.

One of the hardest parts of leaving Michigan will be leaving Mid-Michigan Runners. For
the past few years you have all helped me keep my sanity, and I will really miss Tuesday
night runs. It amazes me how many memories we have made through training runs, the
sweater run, races, vacations, parties, and general get togethers. Thank you, Seth, Simon,
Robert, Lewis, Jenn, and the rest!

Ben, Sam, Joe, Terrence, and Will, you guys are some of my oldest friends. Thank you for
making an eﬀort just about every year for the past 10 years to come see me on my birthday.

iv

I hope we can keep at least a yearly get together going in the future. I want to thank you
guys, and the rest of the Watertown crew, for the much needed breaks throughout graduate
school.

Finally, I want to thank my family. In particular, for the emotional and moral support of
my mother, Sarah Rabin, brother, Jake Butts, and sister, Rachel Butts, that kept me moving
forward throughout my PhD. Thank for being available, especially through the hardest times
in graduate school. And thank you for putting up with my stress and complaining for at least
the past six years! To my girlfriend Fran, thank you for your support throughout graduate
school. You listened to endless rants when deadlines were closing in, but you helped stay
focused on getting things done. Most importantly, thank you for being willing to start a new
job in a diﬀerent state to avoid living across the country from me. I am looking forward to
starting a new chapter in life.

v

TABLE OF CONTENTS

CHAPTER 1:

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 2: DATA-DRIVEN AGENT-BASED MODEL BUILDING FOR
ANIMAL MOVEMENT THROUGH EXPLORATORY DATA
ANALYSIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

6

CHAPTER 3: A GLOBAL SENSITIVITY ANALYSIS OF AN AGENT-BASED

MODEL OF CHRONIC WASTING DISEASE . . . . . . . . . . . . .

33

CHAPTER 4: MATHEMATICAL MODELING OF DISINFORMATION AND

EFFECTIVENESS OF MITIGATION POLICIES . . . . . . . . . . .

40

CHAPTER 5: AN ATTRIBUTE-PRESERVING METHOD TO MINIATURIZE

LARGE SOCIAL NETWORKS . . . . . . . . . . . . . . . . . . . . .

60

CHAPTER 6: STOCHASTIC DIFFERENTIAL EQUATION BASED MODELING

OF DETERRENCE . . . . . . . . . . . . . . . . . . . . . . . . . . .

73

CHAPTER 7: APPLYING REINFORCEMENT LEARNING TO AGENT-BASED

SIMULATIONS OF CONFLICT . . . . . . . . . . . . . . . . . . . .

93

CHAPTER 8: CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

APPENDIX A: SUPPLEMENT TO CHAPTER 2 . . . . . . . . . . . . . . . . . . . 128

APPENDIX B: SUPPLEMENT TO CHAPTER 4 . . . . . . . . . . . . . . . . . . . 136

vi

CHAPTER 1:

INTRODUCTION

1.1 Background on agent-based models
Models are indispensable tools for researchers, oﬀering invaluable insights into various phe-
nomena when real-world experiments are impractical or unfeasible. There are several reasons
that might sway a researcher toward modeling over direct experimentation. For instance, the
sheer cost of certain experiments can be prohibitively high, as exempliﬁed by the expenses
tied to launching payloads into space, which can cost up to tens of thousands of dollars per
kilogram [1]. In such scenarios, modeling acts as a precursor, allowing researchers to assess
the feasibility of their endeavors before actual implementation. Equally, time-critical situa-
tions, such as the planning of evacuation routes amidst looming weather threats, necessitate
immediate action; waiting for real-time data, in such cases, is not only impractical but also
perilous. Beyond these tangible constraints, there also exist ethical barriers. For example,
deliberately releasing pathogens to study the propagation of diseases is not an option. While
these examples are far from exhaustive, they underscore the myriad challenges that models
help researchers navigate.

Just as there are numerous reasons for utilizing models, the realm of modeling is vast.
There are many tools and methods tailored to diverse phenomena. These range from mathe-
matical constructs such as partial diﬀerential equations, to tangible replicas like solar system
models, or even surrogate organisms like rats in drug trials. Among these, computational
models – those that rely on computers to be evaluated – are widely applied across var-
ious ﬁelds, and can signiﬁcantly mitigate the ﬁnancial, temporal, and ethical limitations
associated with conducting experiments. One subtype, the agent-based model, stands out,
particularly for simulating complex, multi-agent systems. This dissertation aims to delve
deeper into this computational paradigm.

What is commonly referred to as agent-based modeling encompasses a trio of model
types: cellular automata (CA), individual-based models (IBM), and agent-based models
(ABM). Of these, CAs, with their simple arrays and iterative rule-based updates, are the
most straightforward. IBMs and ABM, in contrast, encapsulate more complexity. They can
comprise arrays, but their deﬁning feature is the incorporation of agents, entities capable
of mobility and decision-making. Historically, the distinction lay in their focus: IBMs em-
phasized the variability among individuals, whereas ABMs prioritized decision-making and
adaptive behavior. Now, the terms “individual-based” and “agent-based” are often used
interchangeably [2].

1

Fundamentally, agent-based modeling is a simulation technique.

It characterizes the
evolution of a system’s state through a rule-driven process. Here, the “system” encapsulates
the entirety of the model, speciﬁcally both the agents and the environment. A “state”
captures a speciﬁc conﬁguration of these agents and the environment, and “rules” guide the
transitions between these states. The complexity of rules isn’t strictly tied to mathematical
formulations [3, 4]. To illustrate, consider binary CA models such as: John von Neumann’s
self-replicating device [5], James Conway’s game of life [6], or Stephan Wolfram’s Wolfram
models [7]. In these models, their systems consists of single arrays with cells that can have
a value of either 0 or 1. With each iteration, each cell either maintains or switches its based
on neighboring cell values. Even such simple systems can generate many states. For an
array with n cells, there are a total of 2n possible. This complexity is only ampliﬁed in
intricate models. Take, for instance, the forest management ABM described in [8]. This
model contains three types of agents: government, conservationist, and logging-company
agents, each of which has a distinct forestry goal. The model represents the forest has two
two-dimensional array that encodes the availability and monetary value of trees. Each agent
type uses these availability and value maps to construct conceptual maps that encode plans
to fulﬁll their forestry goals. Capturing the vast plethora of potential states in such models
reveals the depth and adaptability of agent-based modeling.

By describing models as transitions between states due to actions being taken, agent-
based modeling employs a bottom-up approach. Such an approach allows for the observa-
tion of emergent collective behaviors among agents in the model. Moreover, this approach
ensures modularity, making complexity additions to the model seamless. As an example, con-
sider modeling an inﬂuenza-like illness. A common modeling approach is to use susceptible-
infected-recovered (SIR) compartmental models [9]; however, agent-based models oﬀer richer
dynamics. Figure 1.1 contrasts an agent-based SIR model with its ordinary diﬀerential equa-
tion (ODE) counterpart. The right panel plots the number of susceptible (green), infected
(red), and recovered (yellow) individuals versus time. The thicker lines are the solutions to
the ODEs and each thinner line is from a single execution of the ABM; in total the ABM
was executed 100 times.

In the ABM version, individual agents move randomly in space. Infected agents, repre-
sented with red points, can transmit the disease to susceptible individuals, represented with
green points, when they come into close contact. Eventually, infected individuals recover
and become immune; these recovered individuals are shown as yellow points. This agent-
based methodology allows for a direct exploration of potential correlations in the model.
For instance, it might provide insights into how the initial spatial distribution of infected
individuals can inﬂuence the rate or extent of disease spread. However, when the goal is to

2

Figure 1.1: On the left is a snapshot of an SIR model implemented as an ABM. Here, red
points represent infected individuals who can infect the susceptible green points. Red points
eventually recover and become immune yellow points. The center shows the mean-ﬁeld
limit of this model, which is the ﬁrst order ODE implementation of an SIR model. On the
right individual runs of the ABM are shown with translucent lines, while bold lines are the
solutions of the ODEs. Notably, while the ODE implementation eﬀectively describes average
behavior, it fails to capture certain dynamics like stochastic extinction, where the disease
ceases to spread before infecting the entire population.

summarize the larger-scale dynamics without getting into the speciﬁcs of individual inter-
actions, we often turn to the mean-ﬁeld approximation. For diseases, this might manifest
as ordinary diﬀerential equations (ODEs) that take into account the average behavior of
agents. In the case of our SIR example, the ODEs are a ﬁrst-order approximation of the
ABM and arise in the mean-ﬁeld limit [10] of the model. This limit assumes that pairwise
correlations among agents factor into a product of averages ((cid:104)SI(cid:105) = (cid:104)S(cid:105) (cid:104)I(cid:105)); however, other
approximations could be used to derive higher order sets of equations [11, 10]. For instance,
a second order approximation would require nine total equations to describe each single and
pairwise correlation of S, I, and R. Although ODE approximations simplify the dynamics,
they can sometimes omit essential details present in the full agent-based model. Such ODE
models condense the complex interactions into parameters, like β for disease transmission,
that average out the diverse individual interactions.

Taking the idea further, to account for agent attributes like age or gender, one would
expand the agent-based model rules to factor in these variables. The diﬀerential equation-
based model would then require a new set of equations and parameters for each added
attribute. This could quickly lead to more complex models like partial diﬀerential equations
(PDEs), which are inherently more challenging to solve. Despite this, the strength of agent-
based models lies in their inherent ﬂexibility to integrate added complexities. However, the
trade-oﬀ is computational eﬃciency; ABMs are generally more computationally demanding
than solving ODE models, especially when one needs multiple model runs, such as in training

3

or uncertainty quantiﬁcation scenarios.

Beyond the cases presented, agent-based modeling has applications in a wide range of
ﬁelds including: ﬁnance [12], optimization [13], human behaviors [14], sociology [15], health-
care [16, 17], and many more [18].

1.2 Areas of improvement for agent-based models
As the ﬁeld of agent-based modeling evolves, it grapples with an inherent tension regarding
the appropriate level of model detail. This tension anchors itself in the age-old principles
of parsimony versus accuracy. Some ABM practitioners posit that models should be as
streamlined as possible, capturing only essential features for clarity, generalizability, and to
prevent overﬁtting. This perspective aligns with Grimm et al.’s discussion the Medawar
zone, which represents the optimal range of model complexity for maximum beneﬁt [19, 20].
Epstein [21] further elucidates this viewpoint, highlighting how even a potentially “wrong”
yet simple model can underpin the foundation of ﬁelds, e.g. Hooke’s Law and the Lotka-
Volterra model. Conversely, Edmonds and Moss [22] introduce the “Keep It Descriptive
Stupid” (KIDS) approach, contrasting it against the traditional “Keep It Simple Stupid”
adage. They argue for comprehensive models, and emphasize simpliﬁcation only when it’s
justiﬁable. Central to this debate is the overarching question of how one strikingly captures
the complexity of reality in computational representations—a challenge that has been echoed
in disciplines beyond ABM, reﬂecting broader epistemological concerns about representation,
comprehension, and forecasting [23, 24].

This competition between parsimony and realism illuminates three facets of agent-based
modeling warranting deeper exploration. 1) Data integration: exploring how data can be
harnessed to determine the necessary complexity of a model and, moreover, once a model’s
structure is established, how data can further guide the determination of its functional form
and parameterization. Rand [25] discusses the development and challenges of agent-based
models in an era of big-data. Additionally, many authors have applied a data-driven approach
to a variety of applications including forecasting diseases [26, 27], family formation [28], and
adoption of solar power [29]. 2) Policy evaluation: delving into how agent-based models can
be eﬀectively applied to real-world challenges, thereby bridging the gap between theoretical
modeling and practical problem-solving. Macal provides an extensive overview of ABM’s
various real-world applications, focusing on both its beneﬁts and the associated challenges
[18]. 3) machine learning incorporation: investigating the fusion of machine learning (ML)
techniques with agent-based modeling to aid in deﬁning the model.

In pursuit of exploring, reﬁning, and pushing the boundaries of agent-based modeling, this
dissertation delves into the previously discussed facets of this paradigm. Each subsequent

4

pair of chapters represents unique projects I undertook aimed at advancing the design and
utility of agent-based models across diverse applications. Chapter 2 incorporates previously
published work [30] that showcases how data can be used to develop an agent-based model
of the movement of deer. Through exploratory data analysis, pivotal features from multiple
datasets were gleaned, culminating in a data-consistent random walk model that underpinned
the agent rules. Chapter 3 builds upon this work incorporating disease and population
dynamics into the model. Using data from multiple sources, the model is parameterized, and
a global sensitivity analysis is performed. This analysis provides methods for quantitatively
exploring variance in the models output to variance in the data used as input to the model.
Chapter 4 pivots to creating agent-based models that can be used to make decisions in the
real world, and is based on previously published work [31]. In particular, this chapter focuses
on strategies to curtail the spread of disinformation in social networks. This chapter presents
an ABM simulating the spread of disinformation and evaluates diﬀerent mitigation strategies,
with the goal of informing users of the model of which strategies are viable. In Chapter 5,
I discuss methods for generating simulated social network. These generated networks are
smaller versions of larger social networks, created such that they preserve attributes of the
original network. The capability to generate networks in this fashion will allow for more
experiments to be run while not compromising the results of the experiments from over-
simpliﬁcations. Transitioning to the intersection of ML and ABM in Chapter 6, we start by
proposing new models of deterrence to use as a baseline comparison to ML-infused ABMs.
Chapter 7, the focuses on the inclusion of reinforcement learning (RL) within ABMs, with
the aim of inferring policies directly from an ABM. We explore deterrence in an agent-based
model of capture the ﬂag where players in the game are controlled by RL agents.

To ensure clarity and ease of navigation, each chapter commences with a succinct sum-
mary, spotlighting the principal contributions and discoveries. Through discussing diverse
applications, from deer movement patterns to the spread of disinformation on social net-
works, this dissertation presents a comprehensive investigation into the future of agent-based
modeling.

5

CHAPTER 2:

DATA-DRIVEN AGENT-BASED MODEL BUILDING FOR

ANIMAL MOVEMENT THROUGH EXPLORATORY DATA

ANALYSIS

2.1 Summary
This chapter delves into the initial facet of ABMs: data integration. The work discussed
in this chapter is based on my previously published work in [30]. The primary goal was to
develop a data-driven ABM to model the movements of deer. To achieve this, two distinct
datasets were analyzed: one of deer relocation data and the other consisting of resource
selection functions. Exploratory data analysis revealed an absence of correlations between
these datasets and a discernible non-parametric trend in the relocation data. To address
this trend, the Fused-lasso machine learning technique was employed, enabling the trend’s
detection and subsequent removal.

My major contribution in this work was the development of a random walk model that
utilized non-Gaussian noise. The random walk model was able to accurately produce sim-
ulated data that was consistent with the data used to develop it, and captured movement
features that are absent in more traditional modeling techniques. Such features are im-
portant when considering the management of diseases and populations of animals, which is
discussed in the following chapter.

Introduction

2.2
Understanding how wildlife move and use the landscape has interested wildlife professionals
since the inception of wildlife ecology and management [32]. Researchers have explored
behaviors such as migration [33, 34, 35], resource use [36, 37, 38], space use [39, 40], and
risk avoidance [41, 42]. Knowledge obtained from such research improves our ability to
manage and conserve wildlife [43]. Quantifying movement behavior is aided by incorporating
observations into mathematical models, such as hidden Markov models [44, 45, 46] aimed at
detecting behavior switching, state-space models [47, 48, 49] used to correct measurement
errors and identify behaviors, and models of random walks on potential surfaces [50, 51] that
attempt to tie movements to resource distributions.

As technology for collecting telemetry data advances, higher-quality data are becoming
available. Such additional data provide insights into complex behavior that can be used
to create more realistic models. One way to bridge the gap between data generation and
analysis is through exploratory data analysis (EDA) [52, 53, 54]. EDA aids the development

6

of more realistic models by enabling features of the data to be identiﬁed and allowing model
assumptions to be checked. EDA methods are broad and vary with the application, but many
approaches, including visualizing data, exploring correlations, and generating statistics, can
be applied to any domain. The importance of checking model assumptions was highlighted in
a series of papers that revealed that 16 out of 17 models claiming to provide evidence for Levy
ﬂights were incorrect [55]. Problems with the analyses included misinterpretation of data
[56], the use of inaccurate ﬁtting methods [57], and the assumption of a heavy tail without
testing alternative hypotheses. EDA can mitigate these problems by allowing assumptions
to be veriﬁed and by exposing model weaknesses.

We provide insight about the importance of EDA and present a speciﬁc example of
applying EDA to model development using movement data and resource-selection functions
(RSFs). The rest of the chapter is summarized in Fig. 2.1 and organized as follows. In
the next section, we give an overview of EDA and present our data and EDA approach.
Our analysis begins with visualizing and identifying features that we later either corrected
for or incorporated into our models. We then test for correlations between the movement
data and RSFs. Next, we present a machine-learning technique to remove trends from
our movement data. The section is concluded by introducing copulas and kernel-density
estimates (KDEs) that are used to construct bivariate distributions from the marginals. In
the following section, we develop a Langevin random-walk model with non-Gaussian noise
that models the movement of a single deer. This model is then discretized to create a
multiple-deer agent-based model (ABM) of deer movement with parameters obtained from
our movement data. We compare our Langevin model to three other random-walk models
to illustrate the importance of EDA and its utility for model development. Finally, we show
the results of ABM simulations of 15 deer in three groups to illustrate how parameters can
be sampled from the data and to illustrate the distributions of areas covered. We discuss
group overlap patterns exhibited in our simulations; studying such group overlap patterns
could lead to a better understanding of disease spread.

2.3 Exploratory Data Analysis
In this chapter, we use EDA to suggest a model’s rules and to provide empirical param-
eters for it. As emphasized by Tukey [52], developing models in this order mitigates the
introduction of biases.

In general, we group EDA into two steps. The ﬁrst step is to explore and adjust data
quality [58, 59]. This step can include data transformations (e.g., logarithmic compression),
imputation, smoothing noisy signals, resampling on more convenient grids and data ﬁtting.
In the second step, patterns are sought in the cleaned data. These patterns may appear as

7

Figure 2.1: Summary of chapter organization. Each green box represents a major step taken
in the corresponding section listed in parentheses.

shapes of distributions, summary statistics, correlations or causal relations. More complex
patterns can be discovered using machine-learning techniques that are capable of ﬁnding
patterns that are diﬃcult for humans to ﬁnd. Central to all of these EDA eﬀorts is visual-
ization. Visualization very clearly reveals patterns [60] such as trends (with line plots) and
correlations (with scatter plots); the functional forms of distributions are readily revealed
and quantiﬁed with, for example, histograms, violin plots and box-and-whisker plots. High-
dimensional data can be viewed using parallel plots [61, 62] and/or using dimensionality-
reduction techniques [63, 64]. All of these analysis and visualization techniques are readily
available in most programming languages, such as in Python’s matplotlib [65], scikit-learn
[66] seaborn [67] and yellowbrick [68] libraries. In what follows, we will employ these EDA
steps to analyze a real dataset. As we will see, some EDA results directly impact the choice of
rules for an ABM, whereas other results provide insight and error detection without directly
impacting model formulation.

2.3.1 Initial Data Analysis
The ﬁrst step in EDA is initial data analysis (IDA), which serves to clean and explore data
without impacting model formation or attempting to answer questions. The goals of IDA
are to remove bad data, identify outliers, expose inconsistencies, perform transformations,
assess the relevance of data and begin visualizations.

To illustrate the IDA process, we begin with a speciﬁc dataset relevant to our goal of
characterizing deer movement patterns. We examine data from GPS-collared deer that were
captured, collared, and monitored according to and with approval by State University of

8

New York College of Environmental Science and Forestry Institutional Animal Care and Use
Protocol no. 2005-1. These data tracked the movements of 71 white-tailed deer inhabiting
central New York, USA from 2008 to 2009 and have been published previously [69, 70].

The data collected on each deer consisted of location data, in the form of GPS locations
and turn angles recorded for each deer in 5-hour intervals, together with the age and sex of
each deer. Locations were recorded for between 8 and 600 days for each deer in this study;
the duration of data collection for each deer varied with the life span of the deer and of the
GPS-collar battery, with how long each collar remained on each deer, and with mechanical
error.

RSFs are also included in the dataset (see Fig. 2.5); RSFs are models that yield values
proportional to the probability that a unit of resource will be used [71]. RSFs can be
constructed using many methods [71], including methods that account for observed animal
movements [51]. The RSFs that we examined were created using a step-selection method to
assess resource selection by the collared deer in this study; these RSFs have been described
previously [69] (link). Six RSFs are included in the dataset, corresponding to resource
selection in three seasons for both males and females.

We began our EDA of this dataset by visualizing the data. Fig. 2.2 depicts a single deer’s
recorded GPS locations in universal transverse mercator (UTM) coordinates in a scatter plot
showing spatial locations and in other plots showing UTM coordinates over time. Note that,
although the collar is attached to a single deer, it reﬂects that deer’s individual behaviors
and its interactions with other deer and the environment. This trajectory shows that the
deer spent time in two regions. We refer to these regions as “basins,” and we refer to jumps
between basins as ”basin hops.” In the time-series data for UTM coordinates, basins can be
seen as relatively constant coordinate values at which the animal remained over a substantial
period of time. A basin hop is indicated by a sudden shift in at least one coordinate. The
number of basins visited varied for each deer in the dataset. Movement within basins and
basin-hopping behavior may reﬂect diﬀerent behavioral states of an animal; for example,
slow motion in a basin may indicate foraging, while basin hops may indicate migration or
dispersal [72].

Positions were recorded with non-uniform time intervals, and there were instances of
missing entries due to GPS collar malfunctions. For later convenience, the dataset was
mapped onto a uniform grid with no missing entries. A uniform grid with Ng grid points
was chosen based on the number of measurements in the dataset. The value at each grid
point was found by linearly interpolating between the closest points in the original dataset
on either side of the chosen grid point.

After interpolating the movement data, we examined the jumps j, calculated as the

9

Figure 2.2: Identifying basins and basin hops. On the left, we show a single deer’s recorded
GPS locations in UTM coordinates. There are two discernible regions, which we refer to
as basins. On the right, the two UTM coordinates are plotted separately vs. the ﬁve-hour
timestep. Basins are indicated by relatively constant locations, and a jump between basins,
which we term a basin hop, is seen as a sudden shift in at least one coordinate.

distances between two adjacent recorded locations. In Fig. 2.3, we show an example of a jump
distribution, which, like the others, closely resembles the zero-mean Laplace distribution

f (j, b) =

1
2b

e−|j|/b.

(2.1)

The parameter b characterizes the scale of the distribution and, for n jumps, has the maxi-
mum likelihood estimate

ˆb =

1
n

n
(cid:88)

i=0

|ji|.

(2.2)

Jump distributions are typically assumed to be Gaussian or power-law distributions [44, 45,
46, 47, 48, 49, 50, 51, 56, 57, 55], resulting in normal random-walk or Levy-ﬂight models,
respectively. A Laplace distribution decays more slowly than a Gaussian distribution, al-
lowing for larger jump sizes; however, Laplace jumps are not as extreme as those of a Levy
ﬂight. Thus, jump sizes allowed by Laplace noise range between those predicted by a normal
random-walk model and those of a Levy-ﬂight model. Comparisons of maximum likelihood
estimates for each type of distribution to the jump data are made in Fig. 2.3, which shows
that a Laplace distribution is superior to a Gaussian, without over-predicting the tail as
would a power law.

We used Akaike’s information criterion (AIC) as a metric to compare Gaussian and
Laplace models for all of our movement data, as shown in Fig. 2.4. For each deer, we treat
UTM Easting and Northing movement data separately, and we compute AIC scores for each
using both the Gaussian and Laplace models. Each blue point corresponds to a single deer’s

10

393000394000395000396000397000UTM Easting4.7374.7384.7394.7404.7414.7424.743UTM Northing1e6394000396000UTM Easting02505007501000125015001750Timesteps4.73754.74004.7425UTM Northing1e6Figure 2.3: Comparison of jump distributions. Shown are the data (grey) and maximum
likelihood ﬁts for power-law (green), Gaussian (blue) and Laplace (orange) distributions.
The data and the Laplace distribution have tail jump sizes in between those of the Gaussian
and power-law models.

interpolated UTM Easting AIC scores, and each orange point to its UTM Northing AIC
scores. The diﬀerence between the AIC for the Gaussian model and the AIC for the Laplace
model for each point is shown on the vertical axis, and the Gaussian AIC for each point is
shown on the horizontal axis. The dashed grey line indicates where the vertical axis value is
zero, i.e., where both AIC measurements are equal. Points lying above the gray line represent
movement data for which a lower AIC is obtained using a Laplace-distribution ﬁt. For all
but two of these points, the Laplace ﬁt had a lower AIC. Both of our proposed ﬁts had a
single parameter, which indicates that for all of our data points lying above the grey line,
the Laplace distribution ﬁt the corresponding movement data more accurately.

After visualizing and interpolating our data, we treated the basin-hopping trends we
observed. Many EDA methods and time-series models assume that data are stationary and
account for trends separately. Before further analysis, the basin-hopping trend needed to be
removed. One method for removing trends from animal-movement data is to use a potential
function [50, 73, 74] whose gradient guides the motion of an animal.

These potentials can be used to construct potential surfaces, which can be linked to the
distribution of resources [51]. The RSFs provided to us could possibly be used to generate
a potential surface, but they were created at a discrete spatial resolution, which resulted in
ﬂat regions with sharp boundaries. These discontinuities make it diﬃcult to approximate

11

100050005001000Jump Size (meters)105104103ProbabilityGaussianLaplacePower LawFigure 2.4: AIC comparison for model selection. Each point corresponds the UTM Easting
or Northing movement data of a single deer. Each blue point corresponds to a single deer’s
interpolated UTM Easting AIC scores, and each orange point to its UTM Northing AIC
scores. The diﬀerence between the AIC for the Gaussian model and the AIC for the Laplace
model for each point is shown on the vertical axis, and the Gaussian AIC for each point is
shown on the horizontal axis. The dashed grey line indicates where the vertical axis value
is zero, i.e., where both AIC measurements are equal. Both models had one parameter;
therefore, points above the gray line correspond to data that are ﬁt more accurately with a
Laplace distribution. All but two of our data points corresponded to movement data that
are ﬁt more accurately with Laplace distributions.

the gradient of an RSF; a smoothing transformation is necessary to accurately incorporate
these RSFs into a model. Moreover, such a discontinuous gradient can mask correlations
between an animal’s positions and an RSF. Using a Gaussian ﬁlter, we created smoothed
maps; Fig. 2.5 shows a comparison of the original and smoothed RSFs. The large regions
of constant values in the original map transition between each other much more gradually,
resulting in a smoother gradient.

2.3.2 Examining Movement Relative to the RSFs
Although our RSFs were not intended to inform animal movement, but rather to inform
seasonal habitat use and selection of resources given availability in the surrounding landscape
[75, 69, 76, 37], we investigated whether correlations between RSF values and deer locations
and preferences for moving to relatively higher RSF values exist coincidentally. The results
for deer locations are shown in Fig. 2.6 and for deer movement are given in Table 2.1.
Note that while there appear to be correlations between the RSFs and deer locations and

12

05000100001500020000250003000035000AIC Gaussian02004006008001000120014001600AIC Gaussian - AIC LaplaceUTM EastingUTM NorthingFigure 2.5: Smoothing resource-selection functions (RSFs). On the left, we show a section
of the original RSF with its sharp boundaries. On the right, we show the same section after
smoothing using a Gaussian ﬁlter. This operation created a smoother function. In these
plots, lighter colors represent higher RSF values, and darker colors represent lower values.
We measured whether animals spent more time at positions with higher RSF values and
whether they tended to move towards positions with higher RSF values and compared the
results for the original and smoothed maps.

movements, they are not strong but are very sex dependent; the reasons for this are unknown.
However, there appear to be no correlations between our RSFs and movement steps, as Table
2.1 shows. Thus, while aspects of the RSFs have potentially useful information, the lack of a
clear picture of the role of this data in describing deer movements led us not to incorporate
this data into our model development.

Change in value

Change in value

using original map male

using smoothed map male

higher
lower
equal

female
27.7% 26.2%
27.2% 26.4%
45.1% 47.4%

higher
lower
equal

female
47.6% 48.3%
47.7% 48.4%
4.6% 3.4%

Table 2.1: Correlations between deer movements and RSFs. This table shows how often
deer moved to locations with higher, lower or equal resource values in the original (left) and
smoothed (right) resource maps. These values are calculated using the six RSFs and all deer
movement data. Each deer uses the map which corresponds to its gender and the time of
year it was tracked. A single deer can use all three seasonal maps for a single gender. The
data reveal a lack of preference for movement relative to these RSFs.

13

Figure 2.6: Exploration of correlations between deer locations and RSFs. The top histogram
shows the probabilities of males and females occupying a location with a given value on the
original set of resource maps, and the bottom histogram shows the same probabilities for the
smoothed set of resource maps. The results are very similar for the two sets of maps and do
not provide evidence that deer locations are correlated with higher values on these resource
maps.

2.3.3 Removing Trends From Non-Stationary Time-Series Data
Because the basin-hopping trend was surprisingly uncorrelated with our RSF values, we were
unable to remove the basin-hopping trend using these RSFs. We then turned to regression
analysis. When a trend has a known functional form (e.g., exponential growth or seasonality),
it can be ﬁt and subtracted from the data. In Fig. 2.2, we see that movement trends in the
time-series data for deer locations are functionally ﬂat regions connected by discontinuities.
The fused-lasso (least absolute shrinkage and selection operator) machine-learning technique
[77] allows us to automatically ﬁnd such a function.

The fused-lasso method aims to ﬁnd a non-parametric function θ that ﬁts uniformly
sampled time-series data y, where y is either of the UTM coordinates with n values (see Fig.

14

0.00.20.40.60.8RSF Value02468Relative ProbabilityOriginalFemaleMale0.00.20.40.60.8RSF Value0246810Relative ProbabilitySmoothedFemaleMale2.2). Because θ can take on any value at any point in its domain, the fused-lasso method
2 = (cid:80)(θi − yi)2.
ﬁrst constrains the ﬁt to be close to the data through minimizing ||θ − y||2
Letting

D =



1 −1
. . .



0

...



0 . . .

0
. . .
. . .
0










. . .

0
...
. . .
0
1 −1

∈ Rn−1×n,

(2.3)

the regularizer λ||Dθ||1 = λ (cid:80) |θi −θi+1| penalizes jumps in the ﬁt θ by pushing the diﬀerence
in adjacent points in θ to 0. Notice that this regularizer has a diﬀerent purpose than it would
in cross-validation. The strength of the penalty is controlled by the hyperparameter λ ≥ 0.
When λ = 0, there is no penalty for jumps, and the resulting ﬁt θ is equal to the data. As
λ → ∞, no jumps are allowed, and θ will be the average of the data. Written mathematically,
the fused-lasso method aims to minimize the relation

min
θ∈Rn

||y − θ||2

2 − λ||Dθ||1.

(2.4)

For a given λ, the quantity in equation 2.4 can be minimized using optimization software;
here, we used CVXPY [78, 79]. Next, we consider how to choose a value for λ.
To choose the optimal λ, we need to deﬁne the critical jump size, deﬁned as

J = ¯j + 3σj,

(2.5)

that diﬀerentiates between a basin hop and random motion. Here, ¯j is the average jump
size, and σj is the standard deviation of the jump sizes in a trajectory. When optimizing λ,
we will compare the jumps in the data, jy, to jumps produced by the ﬁts, jθ. With these
deﬁnitions, we deﬁne another loss function to be optimized. This is done by ﬁnding the
jumps in the data y, and setting to zero the sizes of any jumps that are less than J . This
process results in jumps given by

jy
i =




|yi+1 − yi|,

|yi+1 − yi| ≥ J

.



0,

|yi+1 − yi| < J

(2.6)

Here, the superscript indicates that the jump is calculated from data, and the subscript enu-
merates the total number of jumps. Our objective is to ﬁnd the value of λbest corresponding
to a ﬁt that best reproduces jy
i . The optimization in Equation 2.4 is solved for n values of

15

λ, (λ1, . . . , λn) and the jumps it produces,

jθ
i = |θi+1 − θi|,

are compared to jy

i . We select λk and corresponding θ such that

||jy − jθ||2

2 =

(cid:88)

i

(jy

i − jθ

i )2

(2.7)

(2.8)

is minimized.

Using synthetic data, we illustrate how the fused-lasso approach can be used to detect
basin hops in Fig. 2.7. The left panel illustrates the threshold for basin hops, which is set
at the mean jump size plus three standard deviations. The middle panel shows that the
model, with this choice of threshold, correctly characterizes the largest jumps in the data as
basin hops. Finally, the right panel shows the ﬁt (orange line) to the synthetic data (blue
dots) and shows that this ﬁt successfully captures the functional form of the trend. To test
the robustness of our ﬁtting procedure, we varied the number and sizes of basin hops and
the amount of noise. Fig. 2.8 shows this comparison; in each row, one parameter is varied,
while the other two are held constant. We found that the fused-lasso method reliably ﬁts
basin-hopping in diﬀerent regimes, and we thus used this method to remove basin-hopping
trends from our data.

Figure 2.7: Example of fused-lasso ﬁtting. Simulated data is used to illustrate our fused-
lasso ﬁtting procedure. In the left panel, we show all jump sizes from the simulated data in
blue, and the critical jump size from Equation 2.5 is shown in green. In the middle panel,
we show the critical jumps from data in blue, and those from the fused-lasso ﬁt in orange.
In the right panel, the simulated location data are shown in blue, with the fused-lasso ﬁt
overlaid in orange.

16

0200400Timesteps0.00.20.40.60.8yt+1ytJumps in Time Series0200400Timesteps0.00.20.40.60.8yt+1ytData and Fit Jump Signals0200400Timesteps0.00.20.40.60.81.0y (arb. units)Time Series DatayFigure 2.8: Testing the robustness of the fused-lasso method using well-controlled simulated
data in arbitrary units. We varied characteristics of simulated time-series data, including
the number and sizes of basin hops and the amount of noise, to assess the robustness of our
lasso-ﬁtting procedure. In each row, one parameter is varied, while the other two are held
constant. In the top row, the basin-hop jump size is varied. In the middle row, the maximum
noise value is varied; the noise added in these plots is sampled from a uniform distribution
between zero and the maximum value listed in the title above each panel. In the bottom
row, the number of basin hops is varied. When not varying, the basin-hop jump size is set
to 1, the maximum value for the noise is set to 0.5, and the number of basin hops is set to
4. Our procedure is able to produce realistic ﬁts to the data in many regimes.

2.3.4 Statistical Inference From Stationary Time-Series Data
Beyond visualizing and detrending data, many techniques can be used to mine useful infor-
mation from time-series data [80]. In this section, we introduce two metrics that we used
to analyze time-series data: autocorrelation functions and the mean-squared displacement.
Calculating autocorrelation functions generated insights about repetitive motion in our data,
though such motion was not a feature we had aimed to reproduce in our model. An exam-
ination of mean-squared displacements provided evidence for basins in our data and thus a
rationale for including basins in our models.

Time Autocorrelation Function

The ﬁrst metric we explored was the time autocorrelation function (ACF). ACFs can be
used in combination with Fourier transforms for pattern discovery in time-series data [81] to
reveal insights and test a model’s assumptions [82].

17

02004001.00.50.00.5PositionJump Size: 102004001510PositionJump Size: 1002004001008060PositionJump Size: 500200400200150100PositionJump Size: 10002004001.00.50.0PositionNoise: Max=002004001.00.50.0PositionNoise: Max=0.2502004001.00.50.00.5PositionNoise: Max=0.750200400101PositionNoise: Max=10200400Timestep1.00.50.00.5PositionBasin Hops: 10200400Timestep1.00.50.00.5PositionBasin Hops: 40200400Timestep3210PositionBasin Hops: 90200400Timestep05PositionBasin Hops: 49The ACF R measures how correlated a time series r(t) is with itself over a time lag. For

discrete location measurements, the ACF can be calculated for a lag k = t(cid:48) − t as

R(k) =

(cid:80)N −k

i=1 (ri − ¯r)(ri+k − ¯r)
i=1(ri − ¯r)2

(cid:80)N

,

(2.9)

where ¯r is the average position. Using the Fourier transform

ˆR(ω) =

(cid:90) ∞

−∞

R(k)e−2πikωdk,

(2.10)

the ACF can be decomposed into its component frequencies. These frequencies measure the
time scale of repetitive behavior in the time series. The power | ˆR(ω)| of each frequency ω, can
be calculated to determine the contribution of each frequency to the total ACF. In practice,
we evaluated a discrete Fourier transform using the fast Fourier transform algorithm. For
many deer, we found that the ACF of their positions primarily took the form of exponential
decay, with high- and low-frequency periodic signals superimposed. Fig. 2.9 shows an
example of the ACF and corresponding power spectra, which highlight the high and low
frequencies in the movement data. The high-frequency peak corresponds to daily repetition
in activity, which is likely circadian in origin. The low-frequency peak consists of fewer power
spectral data points and hence is not as statistically signiﬁcant.

The method presented above assumed that the data were uniformly sampled. Alternative
approaches, such as those employing the Lomb-Scargle periodogram [83, 84], allow irregu-
lar data to be treated directly without interpolation. However, because interpolation was
required for other methods used in this work, we implemented the fast Fourier transform
method in our analysis.

Mean-Squared Displacement

The second metric we examined was the mean-squared displacement (MSD). The MSD is
deﬁned as E[|r(t) − r(0)|2]. The MSD measures how far an animal moves on average and
can be used to classify diﬀusion by the power by which the MSD increases versus time.
In Fig. 2.10, we examined the MSD of our stationary (within-basin movement) and non-
stationary (basin-hopping) data. In the non-stationary data, there was no clear relationship
between time and the MSD. In addition, basin hops caused large shifts in the MSD. After
removing the basin-hopping trend, the MSD was constant in time, as expected. A constant
MSD results when an animal remains in one region; for example, the MSD is approximately
constant if an animal remains in a basin.

18

Figure 2.9: Time autocorrelation functions and their Fourier transforms. ACFs of a deer’s
motion in the x, y, and r = (cid:112)x2 + y2 directions are shown in the left column.
In the
right column, we show the corresponding power spectra of these ACFs. The high-frequency
oscillation in the ACFs is seen as a daily peak in the power spectra, and the lower-frequency
oscillation is captured by a wider peak at around one to three months.

19

0.20.00.20.40.60.81.0ACFX ACF05101520253035PowerX Power Spectrum0.00.20.40.60.81.0ACFY ACF051015202530PowerY Power Spectrum020406080100120140160180200220240Days0.00.51.0ACFR ACF110306090120150180360Days0102030PowerR Power SpectrumFigure 2.10: Mean-squared displacements. We show the mean-squared displacements for
our non-stationary data in the left column and for our stationary data in the right column.
In neither case is there an obvious ta dependence; moreover, the stationary data reﬂect a
“non-diﬀusing” deer, as evidenced by the mean-squared displacement being constant over
all but the longest lags.

2.3.5 Recovering Correlations in Position Data
The above EDA explored our data in each coordinate separately. In doing so, we assumed
that the data in each coordinate were independent and imposed an importance on the co-
ordinate system that our data happened to be measured in. For example, by rotating the
coordinates such that a basin hop is perpendicular to one axis, a trend that was observed
in both coordinates before rotation would be observed in only one of the new, rotated
coordinates. Changing coordinates would aﬀect our fused-lasso regression analyses, jump
distributions, and other metrics examined in our EDA. Therefore, it is important to ensure
one coordinate system is used consistently.

One way to remove the assumption of independence between coordinates is through
correlated random walks [85]; however, these models assume that the resulting bivariate
jump distribution is normally distributed, something we did not observe in our data. To
preserve correlations in the data and to avoid imposing a form on the jump distributions
that we pick, we used two methods to estimate bivariate distributions from the data directly:
copulas and KDEs. Copulas allow us to generate a bivariate distribution from the marginal
distributions, which we know well from our stationary data. Jumps in the trends do not
follow an obvious distribution, so we use KDEs to approximate these instead.

20

0.00.20.40.60.81.01.2MSD Meters21e6X Non-Stationary050000100000150000200000MSD Meters2X Stationary0200400600800Day lags0200000400000600000800000MSD Meters2Y Non-Stationary0200400600800Day lags0250005000075000100000125000150000175000MSD Meters2Y StationaryCopulas

Animal locations are often recorded in arbitrary coordinates that are useful for modelers.
There is no reason to assume that motion in these directions should be independent. There-
fore, any distributions stemming from data in these coordinates should be bivariate to ac-
count for correlations. When the marginal distributions are well known, copulas provide a
way to estimate a bivariate distribution while preserving the observed marginals. A cop-
ula C(·) is a function that joins a multivariate distribution function FU V (u, v) to its one-
dimensional marginals FU (u), FV (v) [86] through the relation

FU V (u, v) = C(FU (u), FV (v)).

(2.11)

Using copulas, one can construct a joint distribution of two variables while knowing only the
individual marginal distributions [87]. Thus, we can estimate the bivariate jump distribution
for our data from the one-dimensional distributions we observed. Copulas have been used
in ecological studies to accommodate under-reporting in wildlife-vehicle crash data [88], to
approximate joint space use among animals [89], and to describe distributions of multiple
species of animals [90].

Mappings from the marginals to the joint distribution function are not unique; a C must
be chosen for each application. While there are many options for creating copulas [86], for
convenience, we employed the Gaussian copula by allowing C to be a multivariate Gaussian,
to model jumps in our stationary data (data with trends removed).

In Fig. 2.11, we graphically show an example of the procedure we employed to generate
a bivariate jump distribution using a copula. First, for each individual deer, we identiﬁed
jumps in their separate UTM coordinates. We considered only jumps that happened simul-
taneously in both coordinates, and we treated jumps that happened simultaneously in both
coordinates as single, two-dimensional data points. The resulting two-dimensional distribu-
tion of within-basin jump data was ﬁt with a zero-mean bivariate Gaussian distribution, as
shown in the upper-left panel of Fig. 2.11. The marginals of these data are also shown in
the upper-left panel of Fig. 2.11. This bivariate Gaussian distribution was then sampled for
illustration purposes; the resulting data, together with their associated marginals, are shown
in the upper-right panel of Fig. 2.11. The correlations from our data in the upper-left panel
are preserved here, but the marginals which were Laplace distributions in the upper-left are
Gaussian distributions in the upper-right. Using the probability integral transformation,
the data from the bivariate Gaussian distribution were transformed to have uniform dis-
tributions, while retaining correlations; the resulting uniform, bivariate distribution, shown
in lower-left panel of Fig. 2.11, is the copula of the two jump-coordinate variables of the

21

ﬁtted bivariate distribution shown in the upper-right panel of Fig. 2.11. Finally, the copula
was used together with the observed Laplace-distributed marginals to create a non-Gaussian
bivariate jump distribution with marginal distributions that were constructed from our data.
This bivariate jump distribution is correlated, has marginals that match our data, and can
be sampled eﬃciently.

Figure 2.11: Illustration of a procedure for generating a bivariate probability distribution
from one-dimensional marginal distributions using a copula. In the upper-left panel, jump
data for an individual deer are shown, and the marginals from these data are shown on
the top and right sides of the main ﬁgure in the panel. A bivariate Gaussian was ﬁt to
the marginals data, and this ﬁt is indicated by the orange oval. In the upper-right panel,
points sampled from the the bivariate Gaussian ﬁt are shown, together with their marginals.
Using the probability integral transformation, the bivariate Gaussian data and marginals
were transformed to uniform data and marginals, while retaining correlations; the resulting
uniform, bivariate distribution is the copula and is shown in the lower-left panel. Finally, the
copula was used with the observed Laplace-distributed marginals to create a non-Gaussian
bivariate jump distribution whose marginal distributions are constructed from our real data.

Kernel-Density Estimates

The power of copulas comes from their ability to constrain the marginal distributions of
an approximate bivariate distribution by exploiting known marginal distributions. When

22

known marginal distributions are unavailable, a KDE can be used to estimate a multivariate
distribution ˆf . KDEs have been used in animal ecology studies in methods for estimating
home ranges for animals [91, 92].

To construct a KDE, we consider a probability density kernel K over each observed point
x and create a weighted average of the distances between the current observed point and all
other points in the dataset Xi, through

ˆf (x) =

1
nh2

n
(cid:88)

i=1

K

(cid:26) (x − Xi)
h

(cid:27)

.

(2.12)

Here, h is a smoothing parameter that controls how closely the output resembles the kernel
used. In our application, each x is a two-dimensional point that describes the size of a basin
hop in each coordinate.

We used KDEs to estimate the bivariate jump distribution of the trends removed from
our movement data. First, we identiﬁed the jumps in the trend. In some cases, the jumps
were not aligned in the separate coordinates, and a smoothing procedure was therefore used
to correct for slight misalignments within the data. We approximated the bivariate jump
distribution of the trend using a Gaussian KDE with h = ¯∆r/4 meters, where ¯∆r is the
average distance between points. This value of h was chosen to yield a smooth distribution
function without losing details present in the original dataset.

In summary, throughout our EDA, we explored many aspects of our data to illuminate
features that then guided the selection of the form of our model. Visualizing our data
informed us of the basin-hopping features in our data and demonstrated that jumps in the
data are not normally distributed. We found that correlations between our movement data
and our RSFs are weak, and we have thus excluded RSF data from our model.

When examining the MSDs of our data, we found additional evidence of basins in which
deer locations were constrained. Using these ﬁndings, in the next section, we propose a
model that can reproduce the features of the data revealed in our EDA.

2.4 Modeling
As part of the process of constructing an ABM, we ﬁrst developed a mathematical model of a
single deer’s movement that was consistent with our EDA; this model was then discretized to
form the basis of our ABM. This intermediate step, of generating a model for an individual
deer’s movements before formulating an ABM, allows us to exploit the EDA performed
above.

Our EDA results suggest that a stochastic drift (i.e., basin-hopping), non-Gaussian
random-walk model would capture important features of our data well; hence, we chose

23

to develop a model of this form as our single-animal movement model. Our starting point
was the usual Langevin equation, which is usually written in terms of particle velocities as

dv
dt

= −γv(t) + ψ(t).

(2.13)

Our modiﬁcation of the Langevin model constructs a random process directly from the
data, using a copula to approximate a bivariate distribution of jumps. This results in a
correlated random-walk model with non-Gaussian noise. Moreover, we include basins in
our model through the use of a harmonic potential to constrain deer within a region. To
simulate the basin-hopping trend, basins move according to a Poisson process; this choice
will be discussed in the subsection below.

While our Langevin model is applied to a single animal and misses the eﬀects of hetero-
geneity, these eﬀects can be incorporated into an ABM extension of the Langevin model.
An ABM can model multiple individuals with diﬀerent ages and sexes and can integrate
other data streams, such as knowledge of the locations of bodies of water and impenetrable
objects. Additionally, other features, such as a disease model, are readily incorporated into
an ABM to study disease transmission, without the need to rederive the movement model.
Our goal is for our ABM to serve as a foundational framework into which additional wildlife
models can be incorporated; such a model will facilitate the study of wildlife phenomena
in which movement plays an important, underlying role. In the next two subsections, we
describe the formulation of our Langevin model and then its extension to an ABM.

2.4.1 Single-Animal Langevin Dynamics on a Potential-Energy

Surface

We develop the stationary and non-stationary components of our Langevin model of a single
animal’s movements separately. We will ﬁrst discuss the stationary portion of our Langevin
model. Consider a single tagged animal at a position (x, y) and moving in a two-dimensional
space. The lack of a trend in the MSD data suggests that deer remain localized in a basin
with center (cx, cy) that acts as a potential energy surface, which we model as

U =

(x − cx)2
ax

+

(y − cy)2
ay

,

(2.14)

where ax and ay describe the widths of the basin. The stationary movement within the basin
is modeled by

(cid:34)

(cid:35)

dx
dy

= 2

(cid:34)

γxax
0

0
γyay

(cid:35)−1 (cid:34)

(cid:35)

x − cx
y − cy

dt +

(cid:34) λxx
γx
λyx
γy

(cid:35) (cid:34)

dΨx
dΨy

(cid:35)

,

λxy
γx
λyy
γy

(2.15)

24

which is analogous to a harmonically trapped particle obeying the Langevin equation [93].
The ﬁrst term is the gradient of equation 2.14 and models the eﬀects of the basin, and the
second term contains all of the correlations and describes the random motion of an animal.
Next, we model the animal’s non-stationary basin-hopping. We assumed that basin hops
were independent events that occurred at distinct, random times, and we calculated the
average rate of basin hops from the data. Our assumptions are consistent with a Poisson
process; thus, we used a Poisson process to model basin-hopping behavior.
In a Poisson
process, the inter-arrival times t between k basin hops are described by an exponential
distribution

P (c moves at t) ∼

e−(k/T )t,

(2.16)

k
T

where T is the total number of timesteps at which the movement of an animal is observed,
and k/T is the average rate of basin hops. The cumulative sum of the inter-arrival times
provides the times at which a basin moves instantaneously. The sizes of the basin hops are
sampled from a bivariate distribution constructed from the removed trends using a KDE.

The stationary and non-stationary components of our model are ﬁt separately. In order
to ﬁt the stationary component described by Equation 2.15, we ﬁrst need to remove the
eﬀects of the basin from our stationary data. We do not know the basin’s parameters a
priori, so we choose many possible sets of basin parameters that are then used to subtract
the basin from our stationary data. In practice, we let C1x = 2/λxax and C1y = 2/λyay
and perform a grid search over values for C1x and C1y. Once the basin is removed, we ﬁt
a Gaussian copula, using the method described in Section 2.3.5, to the resulting data to
capture the parameters and random processes in Equation 2.15. To determine the best set
of C1x and C1y, we run our model 10 times for each set of parameters and corresponding
random process ﬁt, and we select the parameters that minimize the least-squares diﬀerences
in the variances and covariances of locations between the simulated location data and real
location data.

To ﬁt the parameters of the non-stationary Poisson process, we set k to the number of
critical jumps observed in a trajectory; then T /k is the ratio of the total number of timesteps
to the number of critical jumps. Each of these parameters, as well as the removed trends,
were stored for each deer. Using the method described in Section 2.3.5, we construct a
bivariate jump distribution from the removed trend with a Gaussian KDE, from which we
sample basin hops.

After ﬁtting the parameters of our Langevin model to all of our deer, we simulated a
trajectory by selecting a set of parameters for a single deer and the corresponding trend and
then running our model. An animal was given an initial position (x(0), y(0)), and the basin
center (cx, cy) was chosen. The entire stationary component of the trajectory was simulated

25

using

(cid:34)

x(t + 1)
y(t + 1)

(cid:35)

(cid:34)

=

x(t) − C1x(x(t) − cx) + Kx
y(t) − C1y(y(t) − cy) + Ky

(cid:35)

.

(2.17)

This equation is a discretized version of equation 2.15. Here, K represents a random number
sampled from the bivariate jump distribution constructed using a Gaussian copula. After
simulating the stationary component, we simulated a trend. This was done by sampling the
times at which basin hops occur, using Equation 2.16.

We then sampled the basin-hop sizes from the KDE constructed from the removed trend
jumps (i.e., the basin hops) and added the simulated trend to the simulated stationary
trajectory. With this model, we are able to simulate the trajectory of a deer with movements
that follow patterns seen in the original dataset.
In the next subsection, we extend this
single-animal model to a multiple-animal ABM.

2.4.2 Multiple-Group Agent-Based Model
We began with a Langevin model that treated a single deer at a time because the data
that we have tracks individual deer. However, in practical applications, we require models
of multiple groups with population diversity. Moreover, many environmental and social
eﬀects are implicitly contained in the harmonic potential used to constrain deer within a
basin. While we would prefer to incorporate these eﬀects explicitly in our models, our
ability to do so is limited without additional data. Nonetheless, to move toward a multiple-
animal simulation, we assumed that a number of non-tracked animals were together with the
tracked animals in a basin and that these non-tracked and tracked animals formed a group
that underwent basin hops together.

Here we provide a simpliﬁed description of our model. For a full description following
description following the Overview, Design concepts and Details (ODD) protocol [94], see
Appendix A. Our ABM consists of group designations (collectives within the population)
and deer. The simulation contains a number of groups, each of which is initialized with a
shared basin center, a number of deer, and a set of parameters and trends trained from a
tracked deer. Based on their parameters, groups determine when they will undergo a basin
hop and how many basin hops they will undergo. At every step of the simulation, the deer in
each group update their positions (x, y) individually using equation 2.17. Each group then
determines whether the current timestep requires a basin hop. If so, the group samples a
jump size from the KDE ﬁt to the trend data and moves the center of the basin along with
its deer according to the sampled jump. Our ABM can be used as a foundation for models
that involve deer movement. For example, a disease or deer-management model could be
added to the ABM to study additional phenomena while allowing for deer movement that is
consistent with real movement data; this will be the subject of future work.

26

2.5 Results
One of our goals was to demonstrate the importance of EDA in developing models. We
compare our Langevin single-deer model (data and basin-hopping; DBH) with three other
models that use diﬀerent assumptions. The ﬁrst two single-deer models are random-walk
models without basin hopping: one with Gaussian noise (Gaussian random-walk; GRW),
and the other with non-Gaussian noise constructed from the location data via a copula
(data and random-walk; DRW). These models ignore trends in the data; such models might
be selected if the data were not visualized and the basin hops were not observed. The
third model incorporates basin hopping and employed Gaussian noise (Gaussian and basin
hopping; GBH). These models are not an exhaustive list of alternatives, but they illuminate
how assumptions made about trends and types of noise can aﬀect the accuracy of movement
models. We additionally present example results from our multiple-deer ABM that display
how groups of deer can come into contact.

The models compared to our Langevin model were ﬁt using methods similar to those
used for our Langevin model. For models without basin hopping, a zero-mean Gaussian
distribution was ﬁt to the bivariate jump distribution constructed from the non-stationary
data. The GRW model sampled this ﬁt directly, while the DRW model passed the same
jumps through a Gaussian copula to simulate a deer’s jumps. Models with basin hops ﬁt a
Gaussian distribution to the bivariate jump distribution of the stationary data. The GBH
model sampled this ﬁt directly, but the DBH model passed the jumps through a Gaussian
copula, so that its jumps were sampled from a Gaussian copula-based distribution. Both the
GBH and DBH models sampled a Gaussian KDE, as described in section 2.4, to simulate a
trend that was used for both models.

We compared the four individual-deer models by training them on a single deer in our
dataset and performing 100 simulations using each model. Both the GRW and DRW models
appear to capture correlations in the location data, as evidenced in Fig. 2.12 by the similarity
of the orientations of the distributions to those of the data. The basin-hopping models
produce very similar distributions in Fig. 2.13 because they share the same basin-hopping
trend in each of the 100 trials, and this trend represents larger jumps than are seen with the
underlying random motion.

Next, we further examined the diﬀerences between models that included jumps sampled
from a Gaussian distribution and those that included jumps sampled from a copula-based
distribution generated using a Gaussian copula. In Fig. 2.14, we show the counts of po-
sitions of a deer moving with Gaussian-distributed jumps within a basin (reds) superim-
posed on the counts of positions moving with copula-based-distributed jumps (blues). Blues
in Fig. 2.14, corresponding to positions moving with copula-based-distributed jumps, are

27

Figure 2.12: A comparison of trajectories for two random-walk models that do not include
basin hops. We compare how often locations were visited using a random-walk model that
uses jumps sampled from a bivariate non-Gaussian distribution constructed using a Gaussian
copula trained on our data, shown in the left panel, with a random-walk model that samples
jumps from a Gaussian distribution ﬁt to our data, shown in the right panel. Each model
was used to generate 100 trajectories, and the densities of the locations visited are shown.

Figure 2.13: A comparison of trajectories for two random-walk models that include basin
hops. We compare a random-walk model using jumps sampled from a bivariate non-Gaussian
distribution constructed using a Gaussian copula trained on our data, shown in the left
panels, with a random-walk model sampling jumps from a Gaussian distribution ﬁt to our
data, shown in the right panels. Both models include the same basin-hopping trend produced
by sampling a KDE trained on our data. Each model was used to generate 100 trajectories,
and the densities of the locations visited are shown. The blue dots show the data on which
the models were trained.

28

spread out over a larger area than reds in the ﬁgure, corresponding to positions moving
with Gaussian-distributed jumps; i.e., the deer movements modeled using a copula-based
distribution covered a larger area than those modeled using Gaussian-distributed jumps.

To further illustrate how these two models diﬀer from each other, in Fig. 2.15, we
show the logarithms of the magnitudes of the diﬀerences between these two sets of position
counts; the position counts generated from Gaussian jumps (reds in Fig. 2.14) are subtracted
from the position counts generated from copula-based jumps (blues in Fig. 2.14). Short-
distance movements are reﬂected in the center of Fig. 2.15, and the green color at the center
of the ﬁgure indicates that far more short-distance movements occur with copula-based-
distributed jumps than with Gaussian-distributed jumps. The yellow points away from the
center of the ﬁgure indicate that occasional larger jumps are also more common among the
copula-based-distributed jumps than among the Gaussian-distributed jumps. Intermediate
distances, shown as red points, were less common among the copula-based-distributed jumps
than among Gaussian-distributed ones.

Figure 2.14: Densities (counts) predicted by both copula-based- and Gaussian-distributed
noise in a basin centered at the origin (cx = cy = 0). Density predictions from our Langevin
model for both Gaussian and copula-based noise, using parameters from a deer in our dataset,
are shown. The deer density for the copula-based model is shown in blues; over that density is
the density from the Gaussian-based model in reds. It is evident that more distant locations
can be reached using a copula-based noise than with a Gaussian noise.

After examining how diﬀerent sets of assumptions aﬀected the movement patterns of our
Langevin-based model, we ran our ABM to explore the behaviors of multiple deer. First, to

29

-1000-50005001000Easting distance from origin (m)-1000-50005001000Northing distance from origin (m)100101102103100101102Figure 2.15: Diﬀerence in deer densities from Fig. 2.14. The log of the magnitude of the
diﬀerences is reported. When using copula-based-distributed movements instead of Gaussian,
there is a large correction for smaller-sized jumps that allow for the occurrence of rare,
very large movements. Further, when using copula-based-distributed jumps, there are fewer
movements to intermediate-range locations compared to Gaussian-distributed jumps.

construct our multiple-deer model, we ﬁt each deer’s data to our Langevin model, creating
distributions of parameters. By sampling these distributions, we created an ABM of three
groups, each containing ﬁve deer. In Fig. 2.16, each group’s approximate density is outlined
with contours in a unique color for each group, with the initial center of the group marked
with an ‘x’ of the same color. The blue group underwent a basin hop into the region occupied
by the green group. This is an example of a behavior that may have been missed using a
model without basin hops. In the context of disease management, this behavior could allow
for the transmission of a disease between two groups.

2.6 Conclusion
We presented an EDA approach as the starting point for model building that resulted in less-
biased models. Our approach was applied to a real dataset consisting of GPS deer locations
and RSFs. The use of an EDA provided a formal testing ground for exploring our data and
testing model hypotheses. For example, we found weak correlations between our movement
data and RSFs, leading us not to incorporate our RSFs in the construction of our models.
It is possible that diﬀerent RSFs, perhaps constructed in a data-driven way to account for
observed animal movements [51], could aid the development of improved, future models.

Our initial EDA identiﬁed important features of our movement data, namely that there

30

-1000-50005001000Easting distance from origin (m)-1000-50005001000Northing distance from origin (m)420246Figure 2.16: Sample multiple-deer agent-based model output. The results of a simulation of
the movements of three deer groups (red, blue and green), each with ﬁve deer, are shown.
The initial location of the center of each group is indicated by a colored ’x’. In this sample,
the blue group underwent a basin hop into the region occupied by the green group.

were disparate geographic regions in which animals spent time, and that animals hopped
between these regions. We referred to these features as basins and basin hops, respectively,
and our models would need to be able to reproduce them. Further exploration of the data
revealed that jumps within these basins do not follow a Gaussian distribution.

In our IDA, we identiﬁed the non-stationary basin-hopping trend that we attempted to
remove from our data with the use of our RSFs. When we found that the movement data and
the RSFs were only weakly correlated, we turned to the fused-lasso machine-learning method
to remove trends from our data. This method successfully and automatically detected trends
in our data and was used to split our data into stationary and non-stationary components.
Beyond our initial data analysis, we examined ACFs and their Fourier transforms, which
highlighted behaviors on multiple time scales. We also calculated MSDs, which provided
additional evidence for the presence of basins that constrained deer movements. We then
used the insights gained from performing our EDA to develop a model of deer movement. We
began with the intermediate step of building a model for a single animal and then extended
this model to a multiple-animal ABM.

31

Our EDA results suggested that a random-walk model that featured a copula-based
jump distribution would capture important features of our movement data. Additionally,
the model would require terms that constrain deer to a basin but allow for the basin to move
intermittently. We created such a model for a single deer and compared it with three other
models that employed assumptions that could be made in the absence of EDA. We found
large diﬀerences in the amount of area covered by deer in each model.

Our single-deer Langevin model was expanded to a multiple-deer ABM that featured
multiple groups of deer. Our ABM was able to reproduce features of deer movement that
our Langevin model could not. In particular, our ABM simulations occasionally revealed
two deer groups moving into the same region. The ability to model such phenomena can be
important for modeling the spread of disease.

In future work, we will extend our movement ABM to model the spread of chronic wasting
disease in white-tailed deer. This disease requires accurate modeling of movement and space
use to understand its transmission and geographic spread.

32

CHAPTER 3:

A GLOBAL SENSITIVITY ANALYSIS OF AN AGENT-BASED

MODEL OF CHRONIC WASTING DISEASE

3.1 Summary
The following chapter builds upon the work discussed in Chapter 2, and is derived from
a paper currently under review, for which I am a coauthor [95]. This work extends the
animal movement model described in [30] by incorporating deer ecology and chronic wasting
disease (CWD) dynamics. The reﬁned model produces both realistic population dynamics of
Midwestern white-tailed deer, and CWD dynamics that are consistent with ﬁeld observations
from Wisconsin. This allows for predictions on the progression of CWD and its impacts on
white-tailed deer populations.

There is a large uncertainty in many of the parameters used in modeling CWD. Because
of this, a comprehensive sensitivity analysis was performed. Such an analysis not only
allows for uncertainty in the results to be reported, but also allows the relative importance
of parameters to be ranked. My main contributions in this project were aiding in the
incorporation of population and disease dynamics, and implementing the sensitivity analysis.

3.2 Background on sensitivity analysis
All models inherently have uncertainties associated with them; such uncertainties may arise
for many reasons, such as unknown initial conditions, variation in parameters, or uncertain-
ties in a model’s structure [96]. Regardless of the form of these uncertainties, sensitivity
analysis is a crucial tool for evaluating the inﬂuence these uncertainties have on the vari-
ability of a model’s outputs [97]. In this chapter, we focus on uncertainties related to model
parameters. Broadly, sensitivity analyses are categorized into local and global types [98].
Local sensitivity analysis explores sensitivity near a nominal point, while global sensitiv-
ity analysis expands this to the entire sample space, including interactions between inputs
[98, 99, 100, 101, 102].

The ﬁrst step of any sensitivity analysis is identifying the quantities of interest, denoted
as Q. These are the metrics that evaluate the model’s performance, and can be direct outputs
from the model or derived from post-processed model data. Q is considered a function of
the p input variables x = [x1, . . . , xp] ∈ Rp. This input vector is a random vector, where
each element is sampled from a corresponding probability distribution.

In local sensitivity analysis, Q is expanded around a nominal point ¯x. While the nominal
point can be any value, commonly it is chosen to be the mean value of the input parameters

33

[96]. Following the notation of [96],

Q(x) = Q(¯x) +

p
(cid:88)

i=1

(xi − ¯xi)

∂Q
∂xi

(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)¯x

+

p
(cid:88)

p
(cid:88)

i=1

j=1

(xi − ¯xi)(xj − ¯xj)
2

∂2Q
∂xi∂xj

(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)¯x

+ . . . ,

(3.1)

where only the ﬁrst few terms of the expansion are written. The second term of the expansion
consists of p derivatives of Q. Each of these derivatives describe the change in Q near ¯x due to
a small change in xi. These derivatives are referred to as ﬁrst-order sensitivities. It is worth
noting, however, that they only measure sensitivity near ¯x, are linear approximations, and
ignore interactions among inputs. To ensure comparability across diﬀerent units, sensitivities
must be scaled. Two scaling approaches are to multiply the ith sensitivity by the mean or
standard deviation of the ith element of ¯x, resulting in the scaled sensitivity coeﬃcient and
sensitivity index, respectively.

A global sensitivity analysis can remedy some of issues of a local one. Speciﬁcally, the full
input space can be explored, the linearity assumptions are eased, and interactions between
inputs can be examined. However, a global sensitivity analysis is often much more involved
to perform compared to a local one [97, 101, 98, 102]. Global sensitivity analysis methods fall
under broad categories, such as regression-based, variance-based, and density-based [103].
In this chapter we will focus on variance-based global sensitivity analysis [104].

Variance-based methods rely on decomposing the variance of the quantities of interest as

var(Q) =

p
(cid:88)

i=1

Di(Q) +

p
(cid:88)

i<j

Dij(Q) + . . . + D12...D(Q),

(3.2)

where Di = var [E(Q|xi)], Dij = var [E(Q|xi, xj)] − Di(Q) − Dj(Q), and so on [101]. From
this decomposition, the sensitivity indices

and

Si =

STi =

Di(Q)
var(Q)
(cid:88)

Sl

l∈I

(3.3)

(3.4)

are calculated, where I is the set of all subsets that include i (e.g. i, ij with j (cid:54)= i, and
so on). Si is referred to as the ﬁrst-order sensitivity coeﬃcient, and STi is the total eﬀect
index [105, 104, 106]. The ﬁrst-order sensitivity coeﬃcient describes the amount of variance
in Q due to the input xi individually, while the total eﬀect index describes the variance
in the output that can be accounted for by xi, including interactions among other inputs.
An overview of performing a variance-based global sensitivity analysis is discussed in [104].

34

This method relies on generating sequences of quasi-random numbers, referred to as Sobol
sampling. Each sequence consists of p+2 points in Rp, where p is the number of inputs to the
model. The ﬁrst two points, A and B, are generated by independently sampling each of the
p parameters. Then, an additional p points are generated by combining components from A
and B, these points are denoted as Ab1, . . . , Abp. A point Abi is generated by replacing the
ith component of A with the ith component of B.

As discussed, it is typically more expensive to perform a global sensitivity analysis com-
pared to a local one, but global analyses are more informative. Mixed methods can alleviate
the issue of computational expense [97] by ranking the importance of inputs using a com-
putationally inexpensive screening method. This ranking is then used to select a subset of
parameters to include in the global sensitivity analysis.

In this chapter, we focus on the Morris method [107] for screening parameters. This
method estimates the ﬁrst-order sensitivities at many random points in the sample space
using a series of trajectories. Each trajectory consists of p + 1 points in Rp. The ﬁrst point
is sampled randomly, then each of the following p points diﬀers from the previous point
in a single element by ∆xi. This process is carried out such that each element is changed
once. Q is calculated at each point in the trajectory (Q(x) and Q(x + ˜xi)), where the
vector ˜xi represents a vector of zeros except for the ith that has a value ∆xi. Using these
measurements, the sensitivity of Q with respect to xi can be estimated using the forward
diﬀerence approximation

dQ
dxi

(x) ≈

Q(x + ˜xi) − Q(x)
∆xi

.

(3.5)

This process is repeated for R trajectories that start at random points in the sample space.
The average, µ, average of the absolute value, µ∗, and standard deviation, σ, of the R
estimates are then used to rank the inputs. The average eﬀect of an input is measured by
µ; however, sign cancellations may lead to misleading conclusions. The total eﬀect of an
input estimate is measured by µ∗ and is unaﬀected by sign cancellations. Higher values of
µ∗ indicate parameters of which Q is more sensitive. Non-linear eﬀects are estimated by σ,
where larger values correspond to more interactions between parameters. Creating a scatter
plot of σ versus µ∗ is a useful tool for choosing important parameters. Parameters that are
close to the origin are less important. Figure 3.1 shows an example with four points. The
parameter x1 would be excluded from further analysis because Q is not sensitive to x1, and
x1 does not interact with other parameters. The parameter x2 does not interact with other
parameters, however Q is sensitive to it. The point x3 interacts with other parameters, but
Q is not very sensitive to it. Lastly, x4 interacts with other parameters, and Q is sensitive
to it. Therefore x2, x3, and x4 would be included in further analysis.

35

Figure 3.1: Example scatter plot of σ versus µ∗. This plot can be used to choose important
features from the Morris method. On the vertical axis, σ approximates how much a pa-
rameter interacts with other parameters. On the horizontal axis, µ∗ estimates how sensitive
Q is to the parameter. In this example, Q is sensitive to the parameters x2 and x4. The
parameters x3 and x4 interact with other parameters. The parameter x1 would be excluded
from further analysis.

3.3 Sensitivity analysis results
In this section, I overview the results of the sensitivity analysis performed on our model of
CWD in Midwestern white-tailed deer. We performed a mixed method sensitivity analysis
that utilized the Morris method to screen 15 parameters, and a variance-based global sen-
sitivity analysis to examine a subset of 11 parameters that the Morris method identiﬁed as
important.

Two modules were added to the movement model, population and disease dynamics.
With the addition of these modules, many parameters were added as well. There is a fair
amount of uncertainty in these parameters. However, they can be tuned such that our model
reproduced realistic short- and long-term population dynamics of Midwestern white-tailed
deer. Additionally, minimum and maximum values for each parameter can be estimated.
Each of these parameters is treated as a random variable.
In an eﬀort to not bias these
parameters, we assume that they were distributed following triangular distributions. Such

36

distributions only require the minimum, maximum, and most likely value of a parameter.
These values are listed in Table 3.1 for each parameter.

name
0) group number
1) direct transmission rate
2) indirect transmission rate
3) prion shedding rate
4) prion half-life
5) disease mortality rate
6) adult male harvest mortality rate
7) adult female harvest mortality rate
8) yearling male harvest mortality rate
9) yearling female harvest mortality rate
10) fawn morality
11) spring immigration/emigration rate
12) fall immigration/emigration rate
13) spring dispersal rate
14) fall dispersal rate

min
50
.0001
.0001
.001
3
.00005
.0002
.0002
.0002
.0002
.0002
0
0
0
0

mode max
251
630
.002
.00052
.002
.00052
1
20
48
120
.002
.00015
.007
.0031
.007
.0013
.007
.003
.007
.00144
.003
.0012
.03
.00267
.03
.00267
.03
.01267
.03
.00433

µ∗
.0074
2321.6
2589.3
0.5952
0.05414
3871.2
698.5
1518.5
654.6
1131.0
1856.5
336.8
468.9
181.8
239.8

σ
.01437
4224.1
4976.8
2.13531
0.15123
10699.0
1401.7
2870.5
1360.9
2579.4
3891.8
724.9
1042.7
329.6
569.4

Table 3.1: Parameter distributions and sensitivity estimates. We list the minimum, mode,
and maximum value for each of the 15 parameters. These values are used to deﬁne a
triangular distribution that are sampled to generate inputs to our model for sensitivity
analyses. The total eﬀect, µ∗, and estimate for non-linear eﬀects σ, are listed in the last two
rows. Using µ∗ and σ as guides, we removed parameters from future analysis (see Fig. 3.2).

Using the distributions for our input parameters, we evaluated 79 trajectories for the
Morris method. Equation 3.5 was calculated 79 times for each input parameter, using a
spacing ∆xi that perturbed the previous point by 10% percent. In total, our model was
run 1264 times to screen the parameters. We approximated µ∗ and σ for each of the input
parameters; they are listed in Table 3.1. In Fig. 3.2, σ and µ∗ for each parameter is shown
as a labeled point. Based on results of the screening, we decided to remove the number of
groups, prion shedding rate, prion half life, and spring dispersal rate from future analysis.
The remaining parameters were used in a variance-based global sensitivity analysis. We
generated 95 Sobol samples from the remaining 11 parameters for a total of 1235 sets of
inputs. In Table 3.2 ,we show the total eﬀect index of our parameters after ﬁve, ten, ﬁfteen,
and twenty years of simulated time.

3.4 Conclusion
Based on our results in Table 3.2, we can see that prevalence of CWD is most sensitive to
harvesting female deer throughout the simulation. Likely, this is due to the fact that when
fawns are born in the simulation, they are placed in a close proximity of their female parent.

37

Figure 3.2: Scatter plot of σ versus µ∗ for the model inputs listed in Table 3.1. The points 0,
4, and 3 are approximately on top of each other. Based on this observation, group number
(0), prion shedding rate (3), prion half life (4), and spring dispersal (13) we removed from our
global sensitivity analysis. We decided to keep fall dispersal rate (14) in our future analysis
as we did not want to full remove dispersal behaviors from the model.

Therefore if that female is infected, or if they are in a contaminated part of the environment,
there is a high probability that the fawn will become infected. These ideas are reinforced by
the low sensitivity to male harvest. However, over time prevalence becomes more sensitive
to harvesting all deer. These observations suggests that culling should be done over longer
timescales, focusing on female deer early on. Early in the simulation, prevalence of CWD is
more sensitive to direct transmission of the disease. As time goes on, prevalence becomes
more sensitive to indirect transmission compared to direct transmission. This is likely due to
the prions needing to accumulate in the environment to have a larger eﬀect. Regardless, the
early importance of direct transmission to late importance of indirect transmission suggest
focusing on culling eﬀorts soon after the detection of the disease. While the results discussed
in this chapter have been speciﬁc to one scenario, our model can be applied to other regions
to model CWD.

38

name
1) direct transmission rate
2) indirect transmission rate
5) disease mortality rate
adult male harvest mortality rate
adult female harvest mortality rate
yearling male harvest mortality rate
yearling female harvest mortality rate
fawn morality
spring immigration/emigration rate
fall immigration/emigration rate
fall dispersal rate

STi(5 yr) STi(10 yr) STi (15 yr) STi(20 yr)

0.073
0.068
0.082
0.056
0.194
0.048
0.102
0.082
0.063
0.126
0.107

0.124
0.065
0.082
0.042
0.186
0.04
0.128
0.091
0.061
0.098
0.084

0.11
0.076
0.073
0.054
0.16
0.056
0.115
0.11
0.062
0.095
0.087

0.084
0.089
0.076
0.078
0.15
0.065
0.123
0.101
0.069
0.093
0.072

Table 3.2: Results from global sensitivity analysis. Here we list the total eﬀect index after
5, 10, 15, and 20 years of simulated time. This index measures the amount of variance in
disease prevalence that can be accounted for by each input, including interactions amount
the inputs. Throughout the simulation, prevalence is most sensitive to adult female harvest
rate and least sensitive to yearling male harvest rate. Early in the simulation, prevalence
is more sensitive to direct transmissions than indirect; however, prevalence becomes more
equally sensitive to these inputs as the simulation goes on.

39

CHAPTER 4:

MATHEMATICAL MODELING OF DISINFORMATION AND

EFFECTIVENESS OF MITIGATION POLICIES

4.1 Summary
The following chapter explores the second facet of ABMs: real-world applications. This
chapter is based on my previously published work in [31]. Here, the focus is on decision-
making processes in the context of disinformation spread on online social networks. A graph-
based ABM was developed and employed to evaluate six strategies for combating the spread
of disinformation. These strategies fall under the broad categories of content moderation,
education, and counter-campaigns, with two strategies from each of these categories being
implemented. This approach involved extensive simulations carried out on thousands of
graph structures, including those modeled after real social networks, to test the eﬃcacy of
these strategies.

My major contributions in this work were implementing the models, and thoroughly
exploring the eﬀectiveness of disinformation combating strategies on various types of graphs.
Care was taken to tie the simulated strategies to real world analogs, ensuring that the ﬁndings
oﬀer valuable insights for practical decision-making in tackling disinformation. This research
not only aids in studying disinformation, but also serves as a guide for creating ABMs that
can aid in decision-making.

Introduction

4.2
The spread of disinformation has brought numerous adverse consequences, such as the
manipulation of the 2016 US presidential election [108, 109], COVID vaccine hesitancy
[110, 111, 112], and the growth of QAnon [113]. Advances in the development of chatbots
are creating new concerns [114, 115, 116, 117]. In response, considerable eﬀorts have been
devoted to detecting [118, 119, 120, 121] and combating [122, 123, 124, 125] disinformation.
Disinformation tracking, bot detection, and credibility scoring tools [126] have been devel-
oped, but the spread of malicious information remains a major challenge. Disinformation has
escalated to a degree that the US Congress is examining intervention policies [127]. These
policies can be divided into two categories: individual-empowering and structure-changing
policies [128].

Individual-empowering policies help individuals to evaluate information they are exposed
Ideally, when
to and include policies such as fact-checking social-media platforms [125].
social-media users interact with disinformation, they would be warned and discouraged from
believing or spreading it further [128]. The reliability of information sources can also be rated

40

[129]. This rating can be performed by experts or users who either rate many articles from
a single source, generating an aggregate score, or rate sources directly. Individuals are more
skeptical of sources with low ratings and are less likely to interact with information from
such sources [129]. In addition to fact-checking and rating sources, one of the most powerful
ways to empower individuals is through education. For example, instructional materials for
teaching critical thinking can be created, such as guides for librarians to teach students to
be aware of fake news [124], or for teachers to teach young students to think critically about
news they come across on social media [122].

In contrast to individual-empowering policies, structure-changing policies prevent individ-
uals from being exposed to certain disinformation entirely. Policies that fall in this category
are primarily implemented by social-network operators who monitor the content on their sites
and remove users or content they deem unacceptable [127]. Additionally, groups can run
counter-campaigns on social media to drown out disinformation with facts [130, 131, 125].

Mathematical models are often used to evaluate, choose and optimize strategies for im-
plementing such policies for combating disinformation, because it is usually impractical to
test such strategies in the real world. Many models can be applied to study the spread of
disinformation [132, 133], including the voter model [134, 135, 136], the Axelrod model [137],
epidemiological models [138, 139, 140, 141, 142], the attraction-repulsion model [143], the
naming-game model [144, 145], and the binary agreement model [146, 147, 148, 149, 150].

In this work, we adopt a policy-driven approach to understanding and combating dis-
information, diverging from the prevalent trend in the literature that often emphasizes the
mathematical or statistical mechanics of models. By both oﬀering the theoretical underpin-
nings of our modeling approach and linking our modeling directly to real-world policies, we
hope to provide insights that can aid in eﬀective policymaking.

Towards this goal, we employ the binary agreement model with committed minorities,
an agent-based model developed by Xie et al. This model was chosen not merely because
it is well studied but because it encapsulates critical properties of disinformation spread.
Furthermore, it highlights a pivotal moment in majority opinion, inﬂuenced by a committed
minority that permits minority rule. In the binary agreement model, a tipping point occurs
when a critical fraction of a network population, pc, which is only approximately one tenth
of the population on complete graphs, advocates a particular opinion strongly [146]. On
heterogeneous graphs, pc is even lower, as the average connectivity of the graph decreases
[151]; such low values of pc make it very challenging to mitigate disinformation. The fact
that this model exhibits a tipping point is important because tipping points have been
observed in real human interactions. For example, in groups assigned to identify an item
in an image, a 25% minority can sway the majority’s answer [148]. Moreover, a goal of

41

disinformation campaigns is often to sway public opinion towards a tipping point. The
ability of a committed minority to overtake majority opinions can have large eﬀects in the
real world, as has been seen with social-media inﬂuence campaigns leading up to the 2016
US presidential election [108, 109], inﬂuence campaigns that have aﬀected societal responses
to the COVID-19 pandemic [110, 111, 112], and the rise of QAnon [113].

The basic binary agreement model of Xie et al.

[146] can be extended to incorporate
more realistic features. Examples of extensions that have been developed include modifying
the propensity of an agent in the mixed state to share one of its opinions [152], varying the
number of interactions needed for an agent to change opinions [152, 153], adding a competing
committed minority [147], expanding the number of possible opinions [154], varying the level
of commitment of the minority [145], and including heterogeneous [155] and dynamic [156]
graphs.

While many studies have examined the eﬀects of a single change to the binary agreement
model on the dynamics of the model, to our knowledge, no study has examined the eﬀects of
incorporating one or more disinformation-management strategies. Here, we explore strategies
that have previously been identiﬁed as potential real-world strategies for combating the
spread of disinformation [127]. We implement individual-empowering and structure-changing
policies – in the form of content-moderation, education, and counter-campaign strategies
– in a modiﬁed version of the binary agreement model, and we examine their ability to
move, smooth, or remove the tipping point exhibited by the model when implemented on
several weighted, heterogeneous networks. We incorporate these strategies into our model
by removing or adding agents or altering agents’ susceptibility to other opinions. Properties
of the tipping point provide us with metrics that can be used to quantify the eﬀectiveness of
various strategies for countering disinformation. We apply our methods to several types of
synthetic networks, including small-world and scale-free networks that capture many features
of social networks. Additionally, we apply our methods to real social-network data, from
Asian users of LastFM [157].

4.3 Methods

4.3.1 Binary Agreement Model and Its Computational Implemen-

tation

We construct agent-based models in which agents are connected by a graph and follow the
binary agreement model proposed by Xie et al. [146] Each agent can hold one of the single
opinions A or B, or the mixed opinion AB, and can be either committed or uncommitted to
that opinion. Pairs of agents update their opinions through one of 12 interactions, using the
update rules given by Xie et al. [146]; for convenience, these update rules are summarized

42

in Table 4.1. We refer to these basic update rules as the opinion update rules.

(speaker

before interaction
opinion
−−−−→ listener)
A A−→ A
A A−→ B
A A−→ AB
B B−→ A
B B−→ B
B B−→ AB
AB A−→ A
AB A−→ B
AB A−→ AB
AB B−→ A
AB B−→ B
AB B−→ AB

after interaction

(speaker – listener)

A − A
A − AB
A − A
B − AB
B − B
B − B
A − A
AB − AB
A − A
AB − AB
B − B
B − B

Table 4.1: Binary agreement model update rules for the three opinions A, B and AB. Here,
our convention is that opinion A is disinformation. Each row of the table shows one possible
interaction.
In the left column, the speaker’s state is listed at the tail of the arrow, the
opinion that speaker shares is above the arrow, and the listener’s state is at the head of the
arrow. The right column shows the outcome of the interaction; the speaker’s state is listed
ﬁrst followed by the listener’s state.

Committed agents are introduced by choosing a fraction, pa, of agents who always hold
the opinion A, regardless of their interactions with others. In our model, this committed
minority spreads disinformation; i.e., the opinion A is disinformation, and the opinion B is
the truth. In general, pa is varied over a wide range to ﬁnd the critical value that deﬁnes the
model’s tipping point; disinformation mitigation strategies aim to move this tipping point
to a more favorable value.

We applied the opinion update rules to agents who were connected on simulated social
networks, which will be discussed in the subsection below. At each time step, for each
agent in the network, one of that agent’s neighbors was selected randomly and uniformly.
For each such pair of agents, either the agent or the neighbor was randomly assigned to be
the speaker, and the other was assigned to be the listener. Probabilistically, the speakers
shared their opinions with the listeners; such interactions determined the new opinions of
both the speakers and the listeners. We show these rules, which govern our computational
implementation of our model, in Fig. 4.1. We refer to these rules as the simulation update
rules.

43

Unless otherwise speciﬁed, all agents not committed to A, i.e., not in the committed
minority, were initialized as uncommitted to B, i.e., as holding the opinion B but not
committed to that opinion. The simulation was updated repeatedly, over many time steps,
until either no agent changed state and none were in a mixed state, or 5,000 time steps had
passed. The fractional densities nA and nB of the nodes was then computed.

We examined the tipping point in our simulations by comparing the fraction of agents
with opinion B (the truth), denoted by nB, with the fraction of agents committed to the
opinion A (disinformation), denoted by pa. Many simulations on each graph type were
performed for each anti-disinformation strategy and the basic binary agreement model; see
the Simulations subsection below for further details. At the end of each simulation, we
recorded the values of nB and pa. We compared these results across simulations; for each
strategy and the basic model, and for each graph type, we plotted nB vs. pa on a single plot.
For the basic binary agreement model without any intervention, a tipping point is evident
on the plot; at the tipping point, nB falls from approximately one to zero once pa ≈ .1, as is
shown by the black dot-dash line in Fig. 4.4 and Fig. 4.5. We compare the locations of the
tipping point, i.e., the value of pa at which nB begins to fall steeply toward zero, for each
strategy to the location of the tipping point in the basic model and to the locations of the
tipping points seen with other strategies.

The goal of our work was to measure the eﬀects of several disinformation-mitigation
strategies on the spread of disinformation using our model. To measure the spread of disin-
formation in our simulations, we examined to what extent each strategy moved, smoothed,
or removed the tipping point exhibited by our model. Concretely, to measure how much a
strategy moved the tipping point relative to another condition, we measured relative changes
in the value of pa at the tipping point between conditions. Smoothing and removal of a tip-
ping point were noted as qualitative changes. These are the outputs of our simulations that
we examine in the Results section below.

4.3.2 Simulated Social Networks
Artiﬁcial social networks were modeled using weighted, directed graphs G = (V, E). Agents
vi ∈ V were represented by the graph’s vertices, and an edge between agents, (vi, vj, wij) ∈
E, represented a connection that allowed vi to interact with vj with probability wij. We
generated several graph types that were initially undirected and unweighted, all of which were
created using NetworkX [158]. The graphs we generated, and their required arguments, are
listed in Table 4.2. We created random and small-world networks by varying the probability
of creating an edge in Erd˝os-R´enyi graphs, and by varying the initial number of edges and the
probability of rewiring an edge in Watts-Strogatz graphs. Scale-free graphs were created by
tuning the number of edges to attach from a new node to existing nodes in Barab´asi-Albert

44

Figure 4.1: Simulation update rules governing computational implementation of the binary
agreement model. At each time step, for each node vi (i.e., each individual) in the graph
(step 1), one of its neighbors vj was selected randomly (step 2). In each pair of nodes, one
was randomly assigned to be the speaker, and the other was assigned to be the listener (step
3). With probability wij (the weight of the edge between the speaker and the listener), the
nodes interacted and updated their opinions (step 4). See Table 4.1 for a table of opinion
update rules.

graphs. Additionally, we examined barbell, lattice, and complete graphs.

Real social networks can be very large. Here, our numerical results are based on relatively
small graphs with 400 nodes, a value chosen to minimize computational cost while retaining
accuracy. This value was chosen by executing the binary agreement model on networks of
progressively larger size and comparing the steady state of the simulations to the mean-
ﬁeld (large-node) limit of the model. Consistent with Xie et al. (see Fig. 1.a) [146], we
found that results obtained with a complete graph with 400 nodes agree well with the mean-
ﬁeld approximation. However, as a ﬁnal study, we performed simulations on the real social
network of Asian users of LastFM [157], which has 7,624 nodes and 27,806 edges; this result
will be discussed in the next section. A graph of the Asian users of LastFM is shown in Fig.
4.2.

The basic topology of G contains the presence or absence of an edge. Each edge was
converted into two weighted, directed edges that formed a loop between the two connected
nodes. These edges can be can be assigned real-valued edge weights that reﬂect the relative
strength of a social interaction. We used a graph’s structure to determine its weights; nodes
were more likely to interact with neighbors that they shared more connections in common
with. Similar to graph transitivity, we deﬁne social transitivity as the ratio of the number of
neighbors shared between the successor and predecessor nodes of an edge to the total number

45

Figure 4.2: Asian users of the LastFM social network. Users are represented as green
nodes, and edges represent a mutual follower relationship between users. This social network
consists of 7,624 nodes and 27,806 edges. The placements of the nodes were determined using
the Yifan Hu algorithm [159].

of neighbors of the predecessor node. We then weighted each edge by its social transitivity.
Mathematically, the weights in the graph are given by

wij =

1 + number of shared neighbors between vi and vj
number of neighbors of vi

.

(4.1)

Here, the additional 1 in the numerator accounts for the fact that vi and vj share a connection.
We used these weights to encode the probability that two nodes will interact in our model.
A listener that shares many neighbors with a speaker has a higher probability of switching
state than one that shares few neighbors with the speaker.

46

graph type
complete graph
watts strogatz graph (small world)
watts strogatz graph (small world)
watts strogatz graph (small world)
watts strogatz graph (random)
watts strogatz graph (random)
watts strogatz graph (random)
erdos renyi graph (random)
erdos renyi graph (random)
erdos renyi graph (random)
barabasi albert graph
barabasi albert graph
barabasi albert graph
grid 2d graph
barbell graph

N
400
400
400
400
400
400
400
400
400
400
400
400
400

k

p m n m1 m2

.02N .5
.12N .5
.25N .5
.02N 1
.12N 1
.25N 1
.02
.12
.25

4
24
50
200

200

199

2

Table 4.2: Graphs that were explored and their parameters. We list the graph types and their
required arguments. For the Watts-Strogatz graphs, the argument k is the initial number of
nearest neighbors that are connected, and p is the probability of rewiring an edge. For the
Erd˝os-R´enyi graphs, p is the probability of creating an edge between a pair of nodes. For the
grid graph, the arguments m and n are the dimensions of the grid. For the barbell graph,
m1 is the number of nodes in each of two fully connected subgraphs that are attached by m2
intermediate nodes. Every graph had 400 nodes, but graphs had diﬀerent numbers of edges.
We also explored a real-world graph of Asian users of LastFM [157], not listed here.

4.3.3 Mitigation Strategies
Using our modiﬁed version of the binary agreement model, we evaluated three mitigation
policies – content moderation, education, and counter-campaigns – that are discussed in the
following subsubsections. In this subsection, we ﬁrst provide an overview of the strategies
we consider for implementing those policies, and the implementation of these strategies in
our model. In the subsubsections below, we discuss details of their implementation in our
model.

We ﬁrst discuss content moderation, in the form of banning some users who spread disin-
formation from social media platforms. We implemented two content-moderation strategies
in our model by removing agents from our graphs in two ways: removing inﬂuential agents
in the committed minority that spreads disinformation, or removing randomly chosen agents
in the committed minority.

Next, we consider two educational strategies. These strategies aim to educate individuals
broadly, and target those who are interacting with disinformation. Broad education consists
of teaching individuals how to identify disinformation to reduce the chances that they will

47

believe disinformation or hold immutable opinions. We refer to this strategy as a skepti-
cism strategy, and we implemented it in our model by reducing committed agents’ level of
commitment to ideas. In contrast, a targeted educational strategy includes fact-checking in-
formation and labeling sources, for example by providing labels on online videos or warnings
on cigarettes, to bias individuals towards the truth. We refer to this strategy as an attentive
strategy, and implemented it in our model by giving all non-committed agents a bias towards
sharing and believing the truth.

Finally, we discuss counter-campaign strategies that counter disinformation with facts.
A counter-campaign strategy can be driven by large agencies; for example, the US Centers
for Disease Control and Prevention provides guidance about sharing vaccine information
[160]. Alternatively, counter-campaigns can be conducted by groups of people who share
information to combat disinformation directly. We modeled counter-campaign strategies by
introducing a second committed minority that is committed to the truth and competes with
the original committed minority that is spreading disinformation.

Using our model, we examined the eﬀects of each of these strategies separately. The
strategies we implemented in our model are shown in Fig. 4.3, and the policies, speciﬁc
strategies, and our implementations of them in our model are given in Table 4.3.

Content Moderation

We examined the strength of content moderation by removing diﬀerent proportions of com-
mitted agents from a graph over an order-of-magnitude range; we removed 0.25%, 0.5%, 1%,
2%, or 2.5% of the total number of nodes in the graph, for values of pa that varied from
0.03 to 0.13. Committed nodes were removed before we executed each simulation examining
content-moderation strategies.

One content-moderation strategy we considered required removing inﬂuential agents in
the committed minority that were spreading disinformation. Thus, we needed to identify
inﬂuential nodes in the committed minority in a graph. We deﬁne inﬂuence using two key
metrics: degree centrality and betweenness centrality. Degree centrality measures the number
of connections an agent has, while betweenness centrality measures how centrally located an
agent is in the graph. Both measures of centrality are normalized to fall between zero and
one. The degree centrality of a node is how many connections that node has divided by the
total number of possible connections, while its betweenness centrality is calculated as how
many shortest paths pass through that node divided by the total number of shortest paths
in the graph. An agent with both high degree centrality and high betweenness centrality has
the potential to spread disinformation directly to many agents and to spread disinformation
to disparate parts of a graph that might not interact without that agent. Using these metrics

48

t
n

i
o
p

g
n

i

p
p

i
t

e
h
t

n
o

t
c
e
ﬀ
e

,
s
e
d
o
n

,
.
e
.
i

,
s
t
n
e
g
a

l
l
a

g
n
i
v
i
g

s
d
r
a
w
o
t

s
a
i
b

a

y
t
i
r
o
n
i
m

d
e
t
t
i

m
m
o
c

e
h
t

n
i

t
o
n

h
t
u
r
t

e
h
t

s
d
r
a
w
o
t

e
z
i
s

d
e
r
i
u
q
e
r

s
e
s
a
e
r
c
n

i

y
l
t
h
g
i
l
s

r
o
f

i

y
t
i
r
o
n
m
d
e
t
t
i

m
m
o
c

f
o

e
u
q
i

n
u

h
t
i
w
s
h
p
a
r
g

n
o

t
p
e
c
x
e

,
y
l
r
a
l
i

m

i
s

m
r
o
f
r
e
p

h
c
a
o
r
p
p
a

l
l
e
b
r
a
b

,
.
g
.
e

,
y
g
o
l
o
p
o
t

d
e
t
e
g
r
a
t

d
n
a
m
o
d
n
a
r

t
n

i
o
p

g
n

i

p
p
i
t

s
e
s
a
e
r
c
n

i

y
l
t
a
e
r
g
m

s
i
c
i
t
p
e
k
s

o
n

o
t

e
l
t
t
i
l

s
a
h

n
o
i
t
n
e
t
t
a

a

f
o

e
z
i
s

d
e
r
i

u
q
e
r

e
h
t

i

y
t
i
r
o
n
m
d
e
t
t
i

m
m
o
c

•

•

•

•

n

i

s
t
n
e
g
a

l
a
i
t
n
e
u
ﬂ
n
i

g
n
i
v
o
m
e
r

,
.
e
.
i

,
y
t
i
r
o
n
i
m
d
e
t
t
i

m
m
o
c

e
h
t

l
a
i
t
n
e
u
ﬂ
n
i

•

s
n
a
b

t
n
e
n
a
m
r
e
p

r
o

y
r
a
r
o
p
m
e
t

s
r
e
s
u

d
e
t
c
e
n
n
o
c

y
l
h
g
i
h

f
o

n
o
i
t
a
m
r
o
f
n
i
s
i
d

g
n
i
d
a
e
r
p
s

n

i

s
e
d
o
n

y
t
i
r
o
n
i
m
-
d
e
t
t
i

m
m
o
c

s
n
a
b

t
n
e
n
a
m
r
e
p

r
o

y
r
a
r
o
p
m
e
t

y
t
i
r
o
n
i
m
-
d
e
t
t
i

m
m
o
c

n
e
s
o
h
c

d
e
t
t
i

m
m
o
c

e
h
t

n
i

s
t
n
e
g
a

y
l
m
o
d
n
a
r

,
.
e
.
i

,
y
t
i
r
o
n
i
m

n
e
s
o
h
c

y
l
m
o
d
n
a
r

g
n
i
v
o
m
e
r

•

h
p
a
r
g

a

n
i

s
e
d
o
n

h
p
a
r
g

a

s
r
e
s
u

n
e
s
o
h
c

y
l
m
o
d
n
a
r

f
o

n
o
i
t
a
m
r
o
f
n
i
s
i
d

g
n
i
d
a
e
r
p
s

e
h
t

n
i

,
s
e
d
o
n

,
.
e
.
i

,
s
t
n
e
g
a

f
o

s
a
e
d
i

o
t

t
n
e
m

t
i

m
m
o
c

f
o

l
e
v
e
l

e
h
t

g
n
i
c
u
d
e
r

y
t
i
r
o
n
i
m
d
e
t
t
i

m
m
o
c

•

•

l
a
c
i
t
i
r
c

g
n
i
h
c
a
e
t

:

m

s
i
c
i
t
p
e
k
s

g
n
i
k
n
i
h
t

r
o
f

s
g
n
i
t
a
r

h
t
i
w
s
l
a
u
d
i
v
i
d
n
i

f
o

g
n
i
k
c
e
h
c
-
t
c
a
f

d
n
a

s
e
c
r
u
o
s

g
n
i
d
i
v
o
r
p

:
n
o
i
t
n
e
t
t
a

n
o
i
t
a
m
r
o
f
n
i

•

•

•

•

n
o
i
t
a
c
u
d
e

49

t
l
u
s
e
r

l
e
d
o
m
n
i

n
o
i
t
a
t
n
e
m
e
l
p
m

i

y
g
e
t
a
r
t
s

y
c
i
l
o
p

n
o
i
t
a
r
e
d
o
m

t
n
e
t
n
o
c

a
p

h
t
u
r
t

e
h
t

o
t

d
e
t
t
i

m
m
o
c

f
o

e
u

l
a
v

g
n

i

p
p
i
t

e
h
t

e
s
a
e
r
c
n

i

,
s
e
d
o
n

,
.
e
.
i

,
s
t
n
e
g
a

f
o

y
t
i
r
o
n
i
m

n
o
i
t
a
m
r
o
f
n
i
s
i
d

r
e
t
n
u
o
c

s
n
g
i
a
p
m
a
c
-
r
e
t
n
u
o
c

r
e
g
r
a
l

•

g
n
i
s
o
p
p
o

n
a

g
n
i
c
u
d
o
r
t
n
i

•

o
t

s
p
u
o
r
g

y
b

s
t
c
a
f

f
o

d
a
e
r
p
s

•

s
n
g
i
a
p
m
a
c
-
r
e
t
n
u
o
c

e
w

s
e
i
c
i
l
o
p

e
h
T

.
l
e
d
o
m

r
u
o

n

i

s
n
o
i
t
a
t
n
e
m
e
l
p
m

i

r
i
e
h
t

d
n
a

n
o
i
t
a
m
r
o
f
n
i
s
i
d

g
n
i
t
a
b
m
o
c

r
o
f

s
e
i
g
e
t
a
r
t
s

d
n
a

s
e
i
c
i
l
o
P

:
3
.
4

e
l
b
a
T

s
e
i
c
i
l
o
p

e
s
e
h
t

g
n

i
t
n
e
m
e
l
p
m

i

r
o
f

s
e
i
g
e
t
a
r
t
s

c
ﬁ
i
c
e
p
s

e
h
t

t
s
i
l

e
w

,
n
m
u
l
o
c

e
l
d
d
i
m
e
h
t

n
I

.
n
m
u
l
o
c

t
s
r
ﬁ

e
h
t

n
i

n
e
v
i
g

e
r
a

d
e
r
e
d
i
s
n
o
c

l
a
n
ﬁ

e
h
t

n
I

.
l
e
d
o
m

r
u
o

n

i

s
e
i
g
e
t
a
r
t
s

e
s
e
h
t

f
o

s
n
o
i
t
a
t
n
e
m
e
l
p
m

i

r
u
o

t
s
i
l

e
w

,
n
m
u
l
o
c

d
r
i
h
t

e
h
t

n
I

.
k
r
o
w
s
i
h
t

n
i

d
e
n
i
m
a
x
e

e
w
t
a
h
t

.
l
e
d
o
m

t
n
e
m
e
e
r
g
a

y
r
a
n
i
b

e
h
t

n
i

t
n
i
o
p

g
n
i
p
p
i
t

e
h
t

n
o

s
a
h

y
c
i
l
o
p

h
c
a
e

t
c
e
ﬀ
e

e
h
t

e
z
i
r
a
m
m
u
s

e
w

,
n
m
u
l
o
c

Figure 4.3: Mitigation strategies examined using our model. Solid arrows represent edges in
the graph, dot-dashed arrows represent a change made to the graph structure, and dashed
arrows represent comments. Blue nodes represent agents with the true opinion, red nodes
represent agents in the committed minority who believe disinformation, and purple nodes
represent agents with variable commitment to an opinion. Two strategies for content mod-
eration, both involving banning users from a network, were implemented in our model by
removing disinformation-spreading, committed-minority nodes from a graph. An education
strategy aimed at increasing skepticism was implemented by reducing the level of com-
mitment to ideas of committed-minority nodes. Another education strategy that brought
attention to disinformation was implemented by biasing all non-committed nodes towards
the truth. We also introduced counter-campaigns into a graph by adding a second committed
minority that was committed to the truth. Each strategy was investigated separately.

and a normalization method, we deﬁne the inﬂuence of a node as:

I(v) =

(cid:18)

1
2

CD(v) + CB(v)

,

(cid:19)

(4.2)

where CD is the degree centrality and CB is the betweenness centrality. This metric allowed
us to identify inﬂuential nodes in the committed minority.
In each simulation in which
inﬂuential committed-minority nodes were removed, the most inﬂuential committed-minority
nodes were removed, up to the desired percentage of all nodes in the graph.

Clearly, if enough nodes were removed, the size of a committed minority would be too
small for it to overtake the majority. Moreover, the size of the committed minority could be
decreased by removing nodes at random as well as by removing inﬂuential nodes. Therefore,

50

we compared the eﬀects of our strategy of removing the most inﬂuential committed-minority
nodes against the eﬀects of a second strategy of removing randomly selected committed-
minority nodes.

Education

In particular, we examined two
Another mitigation policy we examined was education.
strategies based on education: skepticism and attention. We implement a skepticism strategy
in our model by relaxing the strength of commitment, c, held by the committed minority and
endowing committed agents with a probability of changing their opinion after an interaction.
This probability is deﬁned as




p(c) =

2(c − cmin)/(1 − cmin)2

cmin ≤ c ≤ 1,



0

else.

(4.3)

Here, cmin is the lowest level of commitment that a committed individual can have; we
assume cmin = 1/2. Equation 4.3 is a triangular distribution that linearly increases between
zero and four when 1/2 ≤ c ≤ 1. Note that as cmin tends to one, p(c) tends towards a delta
function centered at c = 1. This implies that all agents are fully committed, and thus we
recover the original binary agreement model. In this strategy, a partially committed agent
becomes fully uncommitted if it switches its opinion to the mixed state.

We implemented an attention strategy consisting of fact checking and source rating into
our model by biasing non-committed agents towards the truth.
In our model, the true
opinion is B, and we assume that our attention strategy leads agents to be more likely to
believe B. Therefore, agents in the mixed state are more likely to switch their opinion to
B than to A. Additionally, those who already believe the truth have a preference toward
keeping that opinion. To implement these two eﬀects, we introduce a truth bias, β = .1,
that is the diﬀerence between the probability that an agent with the mixed opinion will
share B in an interaction with a neighbor and the probability that the agent will share
A in an interaction. The truth bias β is also used to increase the probability that an
agent with opinion B will change its opinion after an interaction. Mathematically, we have
P (AB shares B) − P (AB shares A) = β, and P (B switches to AB|neighbor shared A) =
1 − β.

Counter-campaigns

Counter-campaigns were implemented by initializing a group of agents to be committed to
the opinion B, which is the truth, to combat disinformation from those who are committed to

51

A, which is disinformation. To examine the eﬀects of diﬀerent sizes of counter-campaigns, we
considered competing minorities consisting of either 5% or 15% of the total population, and
we initialized these minorities to be committed to B. We assumed that counter-campaigns
would begin in a local region of the graph, similar to an echo-chamber. This eﬀect was
incorporated by initializing the graphs in a non-random way, without changing the topology
of the graph. Using NetworkX’s “spring layout” routine that implements the Fruchterman-
Reingold force-directed algorithm, we assigned a position to each node in the graph that was
contained in the box [−1, −1]×[1, 1]. We then ensured that all nodes that were committed to
A were initialized on the left side of the graph (their ﬁrst coordinate was negative), and that
those that were committed to B were on the right side (their ﬁrst coordinate was positive).
All remaining individuals were initialized as uncommitted to B, i.e., holding the opinion B
but not committed to that opinion.

4.3.4 Simulations
For each of the six anti-disinformation strategies, we conducted 30 simulations on a new
instantiation of each of the 15 graph topologies listed in Table 4.2 and 30 simulations on the
real social network. Each of these 2880 simulations executed our model either 11 (content
moderation) or 15 (education and counter campaigns) times, corresponding to the values
of pa we examined, and had a new initial conﬁguration of committed agents. In total, our
model was executed 118,080 times, and required over 18,000 core hours on a distributed
memory supercomputer. Each simulation was ran until either no agent changed state and
none were in a mixed state, or 5,000 time steps had passed. The results of these simulations
were averaged for each graph topology; this allowed us to examine the variance in outcomes
and outlier results.

4.4 Results
We investigated the eﬀects of three anti-disinformation policies, implemented in the form
of six speciﬁc strategies for combating disinformation, using a binary agreement model we
modiﬁed to incorporate these strategies. We conducted a total of 2880 simulations on 15
diﬀerent graph topologies to explore the eﬀects of each of these strategies, and an additional
30 simulations on a real social network. At the end of each simulation, we recorded the
fraction of agents with opinion B at steady state and the fraction of agents committed to
[147], these values are denoted by
the opinion A. Following the notation in Xie. et al.
nB and pa, respectively. For each anti-disinformation strategy, we plotted nB vs. pa and
examined how the value of pa at the tipping point varied among the strategies employed.
We also examined whether a strategy aﬀected the shape of the tipping point or removed the
tipping point altogether. Below, we discuss our results for each strategy. Because of space

52

limitations, we show results from complete graphs and the real social network in the main
manuscript; the results for all graph types are shown in Appendix B.

4.4.1 Content Moderation
Content moderation was implemented by removing committed nodes from a graph before we
executed our simulations. Committed nodes were removed either based on their inﬂuence
In Fig. 4.4a, we show an example of the eﬀects of content moderation on
or randomly.
a complete graph, and in Fig. 4.5a we show the results on a real social network. Due to
computational limitations, only up to 1% of the committed agents were removed from the
real social network. Results for removing highly inﬂuential committed nodes are shown with
solid lines, and results for removing committed nodes randomly are shown with dashed lines.
Results for the basic binary agreement model in the absence of intervention are indicated
with a black dot-dash line. Diﬀerent line colors indicate diﬀerent percentages of the total
nodes that were removed. Each opaque line shows the average result of 30 simulations,
and each partially transparent line shows the result of a single simulation. Similar plots for
content moderation implemented on all of our graphs are shown in Appendix B.

As expected, as more nodes were removed from a graph, the initial size of the committed
minority needed to be larger to overtake the majority; this feature was shared among all graph
types. For most graphs, removing nodes randomly performed similarly, if not identically, to
removing based on inﬂuence. However, on the real network, removing based on our inﬂuence
metric outperformed removing randomly. As shown in Fig. 4.2, the LastFM social network
has many clusters of nodes. These clusters are interconnected by a few nodes, and our metric
does well at identifying these nodes. For barbell, grid, and the small world graphs with a
low average degree, removing nodes randomly outperformed using our inﬂuence metric. The
unique topology of the grid and barbell network ranks many of the nodes as equally important
(e.g., nodes in the fully connected components of a barbell graph) and continuously ranks
certain nodes as most important (e.g., the center node of the grid graph), while removing
randomly we might ﬁnd a better set of nodes. Once a large number of nodes have been
removed, either method for removing nodes work equally well, except for our real network.
These results suggest using the inﬂuence metric to remove a small number of nodes from
a graph, unless the graph has special topology that may severely constrain how nodes are
selected. If a large number of nodes (relative to the size of the graph) are being removed,
a random approach is better due to the extra computational cost for little to no gain in
performance.

53

4.4.2 Education
We also examined the use of education that aimed to increase individual’s skepticism and
Individuals who are more skeptical of information are less likely to be fully
attention.
committed to an opinion. Whereas attentive individuals are biased towards the truth. Here,
we implemented skepticism by relaxing the commitment of individuals and attention by
introducing a bias towards sharing and holding the opinion B. Using a complete graph in
Fig. 4.4b and a real social network in Fig. 4.5b, we show the eﬀects of skepticism with an
orange dashed line, attention with a green dotted line, a combination of both with a solid
blue line, and no education with a black dot-dashed line. Each opaque line is the average
of 30 runs, and transparent lines are singular runs. We show similar plots for applying
education policies to all graph types in Appendix B.

Attention had a small eﬀect for all graph types. In random graphs, when an attention
strategy was used the size of the committed minority needed to overtake the majority grew
slightly with the average connectivity. Attention had no eﬀect for on the real social network
nor complete, barbell, and grid graphs. Contrary to attention, skepticism required the size
of the committed minority to be between 30% and 35% in order to overtake the majority.
Using both attention and skepticism together produces similar results as just using a strategy
based on skepticism. These results suggest that educating people to be more skeptical can
greatly reduce the dangers of a committed minority.

4.4.3 Counter-Campaigns
Additionally, we examined counter-campaigns that were implemented by introducing an
additional group committed to the opinion B, denoted by pb. For a small counter-campaign,
In Fig. 4.4c, we examine the eﬀects of
we let pb = .05, and pb = 0.15 for a large one.
diﬀerent sized counter-campaigns on a complete graph, and a real social network in Fig.
4.5c. As before, we show the result of no intervention with a black dot-dashed line. We
show the eﬀects of the small counter-campaign with a solid blue line and the large counter-
campaign with a dashed orange line. Each opaque line is the average of 30 simulations, and
transparent ones are individual runs. Again, similar plots for counter-campaigns applied to
all graph types are shown in Appendix B.

As expected, when a larger counter-campaign is used there is a larger portion of individ-
uals with the opinion B at steady-state. The only exceptions are the barbell, grid, and real
social network graphs, where small and large counter-campaigns have the same eﬀect; this
is likely due to their unique topologies. On most graphs when a small counter-campaign is
introduced, the tipping point is located when pa is between .1 − .12, and past the tipping
point nB ≈ 0.05: the size of the counter-campaign. This behavior is not observed on the
grid graph, instead nB appears to decay exponentially to .4.

54

For most graphs, when a large counter-campaign is introduced the tipping point is lo-
cated when pa is between .15 − .2, and past the tipping point nB ≈ .18 − .3. Again,
the exceptions are the barbell, grid and real social network graphs, where large and small
counter-campaigns act similarly and nB decays exponentially to approximately .5. For the
remaining graphs, the ﬁnal proportion of individuals with the true opinion is larger than the
initial counter-campaign, showing the counter-campaign convinced non-committed individ-
uals of their viewpoint. While these results suggest that the larger the counter-campaigns is,
the better, our large counter-campaign was unable to keep the majority of individuals from
holding the disinformation.

4.5 Discussion
In this work, we investigated the spread of disinformation on various weighted, directed,
heterogeneous graphs with committed minorities using the binary agreement model. We
evaluated the eﬀectiveness of three types of policies, namely content moderation, education,
and counter-campaigns. For each policy we implemented two strategies to combat the spread
of disinformation. These strategies were tested on hundreds of graphs based on 15 graph
topologies, and a real social network generated from data of Asian users of LastFM [157].
Over 18,000 core hours were used to explore the eﬀectiveness of the mitigation strategies on
these networks.

Regarding content moderation, we found that removing nodes based on our inﬂuence
metric did not signiﬁcantly outperform removing nodes randomly for most graphs. As more
nodes were removed, the two methods converged towards each other. This is encouraging,
because our simple inﬂuence metric scales as O(|V ||E| + |V 2|), which can quickly become
burdensome to calculate for large graphs. However, on the real social network graph a
more targeted approach outperformed a random one. Therefore, our model suggests that
removing sources of disinformation as they are identiﬁed is a viable method for implementing
content moderation, however a more targeted approach might prevail on graphs with unique
topology.

We also explored education-based policies and found that a strategy that increased peo-
ple’s skepticism had a notably stronger positive impact in our model than a strategy that
focused on people’s attention to disinformation. In fact, our model suggested that strate-
gies such as fact-checking, had little eﬀect, and only showed a small positive eﬀect when
any change occurred. However, it is worth noting that our deﬁnitions for skepticism and
attention are crude due to the simplicity of our model. Therefore, while our model suggests
that strategies that bring attention to disinformation may not be very eﬀective, we caution
against blindly following this advice. We recommend implementing any method that can in-

55

(a) The results of removing the most inﬂuen-
tial nodes are shown with solid lines, and the
results of removing nodes randomly are shown
with dashed lines. Diﬀerent colors indicate dif-
ferent fractions of the total number of nodes that
were removed. Removing nodes randomly per-
forms similarly to removing the most inﬂuential
nodes, as more nodes are removed.

(b) The results of increased skepticism are shown
with an orange dashed line, the results of in-
creased attention are shown with a dotted green
line, and the results of both are shown with a
solid blue line. Attention has little to no eﬀect,
while skepticism greatly reduced the reach of the
committed minority.

(c) The results of a small counter-campaign are
shown with blue lines, and the results of a large
counter-campaign are shown with a dashed or-
ange line. The large counter-campaign was able
to convince more uncommitted individuals of
their opinion, but overall were unable to stop
the minority from taking over.

Figure 4.4: Strategies based on a) content moderation, b) education, and c) counter-
campaigns applied to a complete graph. In each subﬁgure, the fraction of nodes with opinion
B at steady-state, nB, is plotted versus the fraction of nodes committed to A, pa. For refer-
ence, the black dot-dash line shows how nB varies with pA when no nodes are removed, i.e.,
for the basic binary agreement model in the absence of a mitigation strategy. Each opaque
line shows an average of 30 simulations, while the partially transparent lines represent indi-
vidual simulations.

56

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removed0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattention0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.15(a) The results of removing the most inﬂuen-
tial nodes are shown with solid lines, and the
results of removing nodes randomly are shown
with dashed lines. Diﬀerent colors indicate dif-
ferent fractions of the total number of nodes that
were removed. Removing nodes based on inﬂu-
ence outperforms removing nodes randomly.

(b) The results of increased skepticism are shown
with an orange dashed line, the results of in-
creased attention are shown with a dotted green
line, and the results of both are shown with a
solid blue line. Attention has little to no eﬀect,
while skepticism greatly reduced the reach of the
committed minority.

(c) The results of a small counter-campaign are
shown with blue lines, and the results of a large
counter-campaign are shown with a dashed or-
ange line. Both sizes of counter-campaigns re-
sulted in the graph being approximately equally
split in opinions. The larger counter-campaign
resulted in slightly more individuals believing the
truth.

Figure 4.5: Strategies based on a) content moderation, b) education, and c) counter-
In each subﬁgure,
campaigns applied to the Asian users of the LastFM network [157].
the fraction of nodes with opinion B at steady-state, nB, is plotted versus the fraction of
nodes committed to A, pa. For reference, the black dot-dash line shows how nB varies with
pA when no nodes are removed, i.e., for the basic binary agreement model in the absence
of a mitigation strategy. Each opaque line shows the average of 30 simulations, while the
partially transparent lines represent individual simulations.

57

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.8nBnone removed0.25% removed0.5% removed1.0% removed0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattention0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.15crease skepticism, while noting that strategies that focus on immediately biasing individuals
towards the truth may not be as eﬀective as those aimed at more broadly educating people
to be critical thinkers.

After exploring the eﬀectiveness of counter-campaigns, we found that small counter-
campaigns were ineﬀective at combating the spread of disinformation for most graphs. When
employed, they left only the original counter-group holding the true opinion. Contrarily,
larger counter-campaigns were able to sway some uncommitted individuals, suggesting those
committed to the disinformation had less reach. However, in the real social network, grid,
and barbell graphs both large and small counter-campaigns resulted in an approximately
equal split between individuals holding the true or disinformed opinion. We believe this is
due to the topology of these graphs and our initialization of the committed agents. We ini-
tialized counter-campaigns and the original committed minority on opposite sides of graphs
to account for them starting locally on the graph. However, in the real social network, grid,
and barbell graph these two sides are more strongly intraconnected than interconnected,
resembling echo-chambers. Thus each competing group of committed agents overtake their
respective side of the graph, resulting in a fairly equal split of opinions. In contrast, previous
work has explored initializing the competing committed group randomly on graphs [147].

For each strategy we implemented based on content moderation, education, and counter-
campaigns, we found that the fraction of individuals with the true opinion versus the fraction
of individuals committed to the disinformation were quantitatively similar. The vast major-
ity of individuals held the true opinion until a critical number of individuals were committed
to the disinformation. Then there was a tipping point in opinions where the disinformation
group quickly dominated the popular opinion. Past this critical value the number of people
believing each opinion stays constant. The diﬀerent implementations of our polices were able
to quantitatively vary these features. Most notably, education and content moderation in-
creased the required size of the minority needed to sway the majority and counter-campaigns
preventing individuals from being swayed by disinformation. The diﬀerent implementations
of our polices were able to quantitatively vary these features through various mechanisms.
The quantitative changes produced by these strategies are summarized in Table 4.3. Content
moderation and education address two sides of the disinformation issue: the source and the
receiver. Content moderation attempts to curtail the spread at its origin, while education,
especially fostering skepticism, equips individuals to resist misleading narratives. On the
other hand, counter-campaigns introduce a competitive narrative to challenge and diminish
the inﬂuence of the primary disinformation source. Each strategy has its merits, but their
eﬀectiveness can vary based on the structure and nature of the network in question. Most
notably, education and content moderation had stronger eﬀects on the tipping point than

58

content moderation. Speciﬁcally, skepticism moved the tipping point forward, while larger
counter campaigns increased the number of individuals that remained with the true opinion
after the tipping point. As the ability of generative models to produce convincing material
grows they might be used to generate disinformation that is diﬃcult to detect [117]. In such
scenarios, educating individuals will become evermore important.

Based on the results of our simulations, education-based policies that increase skepticism
and counter-campaign policies were the most promising. Examples of real world policies that
increase skepticism could include creating media literacy programs. Such programs teach
individuals to recognize trustworthy sources and discern fact from opinion. Additionally,
incorporating critical thinking curriculum into the education curriculum at an early age can
build a foundation of skepticism. Another potential strategy would be creating feedback
loops on social media platform that notify users when they engage with potential disinfor-
mation. Example of real world strategies that focus on counter-campaigns could range from
launching corrective advertising campaigns that directly address and correct false narratives,
to establishing rapid response teams. These teams would monitor social media for emerging
disinformation so that they could quickly launch counter-eﬀorts.

In summary, our goal in this work was to create a model that captures realistic mitiga-
tion strategies within the conﬁnes of a binary agreement model so that our result may be
applicable to both modelers and policymakers. Additionally, there is a vast number of other
policies, strategy implementations, and combinations of the two that could be explored. This
leaves many avenues for future studies to extend our work. One promising direction is to in-
corporate a more realistic model of the spread of opinions, such as the Attraction-Repulsion
Model [143]. This model allows for multiple, continuous-valued opinions to be held by each
individual. This more nuanced representation of opinions allow individuals to become more
or less closely aligned with their neighbors. Additionally, we could enhance the realism of
the model by incorporating multiple online and oﬄine social networks, dynamic networks,
and news sources. Beyond improving the model, future studies can explore policies we did
not test, optimal parameters for existing policies, or learning entirely new policies.

59

CHAPTER 5:

AN ATTRIBUTE-PRESERVING METHOD TO MINIATURIZE

LARGE SOCIAL NETWORKS

5.1 Summary
This chapter builds upon the work discussed in Chapter 4. Speciﬁcally addressing the
diﬃculty in running simulations on real social networks, as well as the non-representative
In this chapter, a method for miniaturizing social
nature of simulated social networks.
networks is presented; though, this method can be applied to any type of network. There
are many ways to generate graphs, such as Erd˝os-R´enyi and Watts-Strogatz models; however,
this proposed method diﬀers in that it creates graphs that have multiple similar metrics to
a given graph but are reduced in size. Many more simulations can be evaluated on these
miniaturized graphs compared to large network, without sacriﬁcing accuracy.

My main contribution in this chapter is the implementation of a parallelized version of
the parallel tempering algorithm. This algorithm generates miniaturized versions of social
networks, with the constraint that the error in a set of metrics is minimized.
I present
an example of miniaturizing a real social network while maintaining the density, degree
assortativity coeﬃcient, and clustering coeﬃcient. This example provides all the necessary
steps to initialize and execute the algorithm.

Introduction

5.2
Networks serve as a fundamental tool for modeling systems composed of interconnected
entities, oﬀering insights across various disciplines [161]. Biologists, for example, utilize
networks to describe genetic dependencies impacting protein expression [162, 163, 164]. In
the realm of social sciences, networks represent social interactions, aiding in understanding
how community structures inﬂuence behavioral dynamics [165, 166, 167, 168]. Similarly,
networks are fundamental in modeling contagion-like phenomena such as the spread of dis-
eases [169, 170, 171], ideas [29, 172], and misinformation [31, 173]. Networks also shed light
on the intricacies of broader social-ecological systems, guiding scientists, stakeholders, and
policymakers. The proliferation of large network datasets, increasingly available through
public data repositories [174, 175, 176], open new avenues for simulating real-world systems
on realistic graphs. These datasets present notable opportunities to simulate the dynamics
of many real-world systems on realistic graphs.

However, the large size and complexity of existing network datasets pose signiﬁcant scal-
ability challenges, particularly for researchers relying on graphs as a tool for modeling and
analysis. As networks grow in size and complexity, running simulations incur prohibitive

60

computational costs, often restricting our ability to accurately model associated real-world
systems [177]. While new scalable algorithms have been proposed to mitigate these issues
[178, 179], they are often developed for speciﬁc modeling tasks and may not be generally ap-
plicable outside of that task. Additionally, computing solutions that have sought to capital-
ize on the increased capabilities of distributed memory computing clusters require expertise
in graph-based parallel computing paradigms, imposing additional constraints on network
science practitioners [180, 181].

One approach to addressing the scalability challenges of large networks is the generation
of structurally similar yet miniaturized synthetic graphs, enabling more manageable analysis
and simulations. Existing graph generation methods excel at matching a narrow set of
characteristics of the graphs that they represent at the expense of control over a broader
network metrics. For instance, Erd˝os-R´enyi [182] and Barab´asi-Albert [183] models are adept
at replicating speciﬁc graph characteristics, such as degree distributions, but they fall short
in controlling broader network features, such as community structure. Our proposed graph
generator is designed to provide more generalized control over a range of properties, aligning
closely with the structural intricacies of large-scale real-world networks. This approach
not only facilitates eﬃcient model testing but also ensures that the models retain greater
relevance and applicability to real-world scenarios.

We propose a graph generator that employs the optimization technique parallel tempering
[184, 185]. This algorithm consists of multiple instances of the Metropolis Monte Carlo
algorithm, called replicas, each at a unique temperature level. Periodic swaps of temperatures
between these replicas allow the algorithm to overcome shortcomings traditionally associated
with optimization, such as convergence to sub-optimal solutions. The process begins with a
seed network – a smaller, preliminary version of the network – and a set of graph metrics
to match from a larger network. The algorithm then iteratively reﬁnes this seed network,
ensuring that it evolves towards resembling the target graph in terms of the chosen metrics.
This method oﬀers an eﬀective means of generating miniaturized networks that retain the
structural essence of their larger counterparts.

5.3 Methods
Our objective is to generate a miniaturized undirected graph G = (V, E) with vertices v ∈ V ,
and edges (vi, vj) ∈ E, such that it closely resembles a set of graph properties of a given target
graph Gtarget = (Vtarget, Etarget). Importantly, we aim to maintain structural and statistical
similarity while ensuring |V | < |Vtarget| and |E| < |Etarget|, thereby achieving a reduction
in both the number of nodes and edges compared to Gtarget. This reduction is crucial for
simplifying the computational complexity while preserving the essential characteristics of the

61

original network. We employ the Monte Carlo optimization algorithm, parallel tempering,
known for its eﬃcacy in sampling complex optimization landscapes by avoiding local optima
to achieve this goal. The following subsections detail the algorithm and how it is applied to
optimize miniaturized graphs.

5.3.1 Metropolis Monte Carlo
As discussed, our objective is to create a graph that minimizes a sum of metric errors, χ,
for various graph properties like assortativity, density, and clustering coeﬃcient. The error
function is deﬁned as:

(cid:88)

χ =

χi(G, Gtarget),

(5.1)

i

where each χi represents the error in a speciﬁc graph property, G is the graph under con-
sideration, and Gtarget is the target graph with desired properties. The Metropolis Monte
Carlo (MC) algorithm forms the foundation of our approach to minimizing χ. The MC al-
gorithm on its own can be used to explore the space of possible graph conﬁgurations, and is
particularly useful for optimizing graph properties due to its ability to sample from complex
probability distributions.

The algorithm is deﬁned as follows.

Initially, a seed graph G0 is initialized.

In each
iteration of the algorithm, a uniformly random number of modiﬁcations between 1 and
10, such as adding, removing, or rewiring edges, are proposed to the current graph G to
obtain a new graph G(cid:48). In this implementation, moving an edge has a 50% chance of being
selected while adding and removing each have a 25% chance of being selected; however these
probabilities can be adjusted. The change in error, ∆χ = χ(G(cid:48))−χ(G), is calculated, and the
proposed change is accepted with a probability p = min (cid:0)1, e−β∆χ(cid:1), where β represents the
inverse temperature. The process is repeated with a new proposed change until a speciﬁed
number of iterations is reached. While the algorithm converges to an optimal solution given
enough iterations and a suitable choice of β, setting a value for β itself is challenging. When
β is low, many conﬁgurations are explored without converging towards a speciﬁc solution.
Conversely, when β is high, the algorithm will quickly converge to the closest minimum of
χ, regardless if it is optimal. This balancing of β is known as the exploration-exploitation
trade oﬀ. Practitioners suggest tuning β such that the acceptance probability is about 23%
[186, 187] to balance the trade oﬀ.

This balancing can be achieved by assuming

χi =

(cid:88)

i

ωi|mi(G) − mi(Gtarget)|,

(5.2)

where each mi ∈ [0, 1] is a min-max scaled metric that is weighed by ωi. This scaling ensures

62

that metrics are on a comparable scale, preventing any one metric from disproportionately
inﬂuencing the optimization process and makes diﬀerences in weights more directly compa-
rable. The weights ωi are set such that each metric is equally weighted based on how much a
they change on average when a graph is perturbed. This balancing is achieved by executing
a MC simulation with β = 0 and ωi = 1 for Test iterations. After each iteration j, the value
|mi(G) − mi(Gtarget)|j is recorded for each metric. Based on these values, the weights are
adjusted as:

ωi =

(cid:80)

Test
j |mi(G) − mi(Gtarget)|j

,

(5.3)

ensuring a balanced consideration of all metrics. Subsequently, another MC simulation is
executed with these determined weights for another Test iterations. After each iteration, the
absolute change in error (|∆χ|)j = (|χ(G)−χ(G(cid:48))|)j is calculated. From these measurements,
the optimal β is deﬁned as

¯β = −

Test ln(.23)
(cid:80)
j(|∆χ|)j

,

(5.4)

which inversely relates to the acceptance probability. With this optimal value for β, we move
to extending the MC algorithm in the next section.

5.3.2 Parallel Tempering
Even with an optimal β to balance exploration-exploitation trade oﬀ, a single instantiation
of the MC algorithm may still struggle to ﬁnd the global optimal solution in complex er-
ror landscapes. Parallel tempering addresses this challenge by running multiple Metropolis
MC simulations (replicas) simultaneously, each assigned a diﬀerent β. This multi-replica
approach allows hotter (lower β) replicas to explore the surface broadly, while cooler (higher
β) replicas focus on exploiting promising regions. Periodic exchanges of β values between
adjacent replicas prevents replicas from getting trapped in local optima, while maintaining
their ability to converge towards optimal solutions.

The parallel tempering algorithm is deﬁned as follows. Initially, N replicas of the MC
algorithm are initialized, each with a unique βi and initial graph. Each replica independently
runs the MC model as described in Sec. 5.3.1. However, at every Tswap iterations, adjacent
replicas propose to swap their β values. This swap is accepted with probability αswap =
min (cid:0)1, e(βi−βj )(χ(Gj )−χ(Gi))(cid:1), facilitating dynamic adjustment to the exploration process.

We implement a parallelized version of the parallel tempering algorithm using OpenMPI.
In this implementation, each replica runs independently on a separate MPI rank, except
during the swap proposal phase. This approach allows for eﬃcient execution of multiple
replicas with minimal additional computational cost. Communication among the ranks in-
volves a structured multi-step process. We assume an even number of replicas, with ranks

63

sequentially numbered starting from 0. In the ﬁrst step, even-numbered ranks send their
inverse temperature βi and current error χ(Gi) to the next odd-numbered rank i + 1. The
rank i + 1 evaluates αswap to determine if a swap should occur. If a swap is to occur, rank
i + 1 sends rank i its own βi+1 and then both ranks exchange there own value of β for their
neighbor’s value; if not, no exchange happens. In the subsequent step, even-numbered ranks
greater than zero send their βi to the preceding rank i − 1 that then decide whether to swap
β values in a similar manner. A visual schematic of this communication process is illustrated
in Fig. 5.1, providing a clearer understanding of the algorithm’s structure.

Figure 5.1: Schematic of communication between ranks. In the ﬁrst step, even-numbered
ranks send their current βi and χ(Gi) forward. If the corresponding odd rank determines
a switch should occur, it sends βi+1 backward and updates its own βi+1 to βi, otherwise it
sends βi backward. Next, even-numbered ranks greater than 0 send their current βi and
χ(Gi) backward and a similar process is repeated.

In addition to the optimal ¯β calculated in Sec. 5.3.1, we incorporate a range of β values to
achieve diﬀerent average acceptance probabilities. We introduce both higher (more exploita-
tive) and lower (more exploratory) β values than ¯β. Here, we set acceptance probabilities
for adjacent replicas to vary by a factor of two (e.g. ¯β/2, ¯β, 2 ¯β). This strategy is designed
to ensure a broad exploration of the solution space, while still eﬀectively converging to op-

64

timal solutions. Employing the Metropolis Monte Carlo and parallel tempering algorithms
as outlined, we shift to applying these techniques to the miniaturization of a social network,
as detailed in the next section.

5.4 Results
In this section, we present the results of applying our graph miniaturization methods to the
Hamsterster online social network [188, 174]. This original network consists of 2426 nodes and
16,630 edges, representing individual users and their friendship connections, respectively (see
Fig. 5.4). We miniaturized this network to 600 nodes, achieving an approximately fourfold
reduction (see Fig. 5.5). The primary goal of this miniaturization was to closely match the
original network’s density, degree assortativity coeﬃcient, and clustering coeﬃcient.

For the miniaturization process, our seed graph was an Erd˝os-R´enyi graph with 600 nodes
and edge inclusion probability of .01. To determine the weights for the metrics, we conducted
10 independent MC simulations, each running for 100 iterations. These simulations were used
to evaluate Eq. 5.3 for each metric. The computed average weights from these 10 simulations
are listed in Table 5.1. Using these weights, we performed another set of 30 independent
MC simulations, each running for 100 iterations, to calculate ¯β = .78873 per Eq. 5.4.

metric
density
degree assortativity coeﬃcient
clustering coeﬃcient

weight
166081
1146
8709

Table 5.1: Weights (ω) for metrics. These weights were determined by averaging the results
from 10 independent MC simulations and evaluating Eq. 5.3. Using these values will weight
each metric equally.

Using the calculated weights and ¯β, we executed the MC algorithm for 10,000 iterations.
The outcomes of these simulations are shown in Fig. 5.2, where panel (a) illustrates the
trajectory of χ, and panels (b), (c), and (d) display the absolute percent error in the degree
assortativity coeﬃcient, density, and clustering coeﬃcient, respectively. In these ﬁgures, in-
dividual results from the 10 simulations are shown in grey lines, while the black line is the
average trend. Notably, χ decreases smoothly across the simulations, indicating an overall
successful miniaturization, with most metrics trending towards minimal error. However,
matching the clustering coeﬃcient proved to be more challenging, as evidenced by its persis-
tent error rate of over 40% even after 10,000 iterations, in contrast to the near-zero error in
the other metrics. The algorithm more readily captured the target network’s density, show-
ing a smoother and more steady convergence towards zero error, compared to the degree
assortativity coeﬃcient, which also showed improvement but with more variation.

65

(a) χ

(b) degree assortativity coeﬃcient

(c) density

(d) clustering coeﬃcient

Figure 5.2: Loss (χ) and target metrics versus iterations for 30 MC runs. The grey lines are
the results from individual runs, and the black line is the average of the result. Overall, χ is
smoothly decreasing for each of the runs. The degree assortativity coeﬃcient and density are
being matched well on average, though the degree assortativity coeﬃcient has much more
variance. While the clustering coeﬃcient is decreasing, it struggles to get below 50% error.
Even with the challenge of matching clustering, a single MC run is trending towards the
target metrics.

Following from this example we turn to the parallel tempering algorithm. Here we chose
to use 6 replicas. This includes ¯β, three lower and two higher betas. Here the set of βi =
{1.90472, 1.53273, 1.16073, 0.78873, 0.41674, 0.04474}. In Fig. 5.3, we show the trajectory of
χ of each replica and their values for each metric, in addition to the β value each replica
has throughout the simulation. With the exception of the trajectory with the lowest β, χ
is trending downward for all replicas. As with the single replica run, the degree assorativity
coeﬃcient and density are fairly easy to match though with more variance in the degree
assortativity coeﬃcient. Given enough iterations, replicas with higher β values correctly
capture the clustering coeﬃcient as well. Exchanges of β between replicas mostly happen
during isolated parts of the simulation and mostly occur between the central values. The

66

steady decrease in error among the metrics, coupled with the the observation that replicas
with higher β values perform better suggest that the error surface has a dominating minimum
that is fairly easy to ﬁnd. This is further illuminated by the replicas with low β values
struggling to ﬁnd a region where the metrics signiﬁcantly deviate enough to outperform
other replicas.

To highlight the similarity in the miniaturized and target networks structures, we display
the target graph and best miniaturized graph in Fig. 5.4 and Fig. 5.5, respectively. Both of
these graphs were drawn using the force atlas algorithm in Gephi [189]. The target network
consists of 2426 nodes and 16,630 edges, while the miniaturized one consists of 600 nodes
and 1002 edges. Many of the features displayed in these ﬁgures are driven by the method
used to draw them; however, both appear similar in that they have a centralized connected
component surrounded by smaller sub-graphs. In both graphs, the central component also
appears to have long strings of connected nodes.

To further examine how similar the generated graph’s metrics are to the target graph’s
metrics, we list the value of the metrics after 20,000 iterations, as well as their initial values
and the target values in Table 5.2. In agreement with Fig. 5.2, we see that all of the metrics
converge to values near their targets.

metric
degree assortativity coeﬃcient (scaled)

initial
0.46803
degree assortativity coeﬃcient (unscaled) −0.06393
0.01002
0.01003

density
clustering coeﬃcient

ﬁnal
0.52333
0.04667
0.00557
0.53757

target
0.52370
0.04740
0.00565
0.53753

Table 5.2: The initial and ﬁnal metrics of a miniaturized network versus the target metrics.
As displayed in Fig 5.2, each of the metrics eventually converge to their target values.

67

(a) χ

(b) degree assortativity coeﬃcient

(c) density

(d) clustering coeﬃcient

(e) β

Figure 5.3: Loss (χ), target metrics, and β versus iterations for a parallel tempering run,
with each color representing a replica. As with the single replica run, the degree assortativity
coeﬃcient and density are matched quickly, with more variation in the degree assortativity
coeﬃcient. Given enough iterations the cluster coeﬃcient is also matched while preserving
the other metrics. The replica with the smallest β explored many conﬁgurations but did not
share information outside of one exchange.

68

Figure 5.4: Hamsterster social network. Blue nodes represent users of the online social
network, and black edges represent a friendship between users. This graph has 2426 nodes
and 16,630 edges. The attributes we aim to match for this graph are listed in Table 5.2.

69

Figure 5.5: Example of a miniaturized graph. This graph is a miniaturized version of the
Hamsterster social network. Similar to Fig. 5.4, the blue nodes represent users of the social
network and the black edge represent friendships between users.

70

5.5 Discussion
In this work we proposed a method for miniaturizing a target network such that the resulting
network has attributes that resemble the target one. To achieve this, we implemented a ver-
sion of the parallel tempering algorithm. We present a method for selecting the optimal β for
a single replica, and then produce β values for additional replicas from the optimal one. Our
implementation of the parallel tempering algorithm is parallelized in such a way that allows
for additional replicas to be added with minimal computational cost. Using our algorithm,
we presented an example of successfully miniaturizing a graph by approximately fourfold,
while preserving the degree assortativity coeﬃcient, density, and clustering coeﬃcient.

In our example, we were careful to chose a target network that was large enough to be
miniaturized by a fair amount but not so large we could not display it. Another advantage
to our target network, was having the ability to exactly calculate our metrics of interest in a
short time. Nonetheless, our methods can be applied to any network, so long as the target
metrics can be calculated.

In our exploration of miniaturizing graphs a single replica performed fairly well. We
tuned β by running preliminary simulations and adjusting β until it accepted around 23%
of the suggested changes. While mildly tedious, it produced decent results by matching
two of the three metrics and tending towards matching the last. By shifting to parallel
tempering we were able to achieve much better results. However, the main beneﬁt of gaining
exploration did not aid our search. While there were exchanges between the replicas, they
were not among the best performing replicas. We believe this is mostly driven by the metrics
we chose to match causing the error surface to not be very complex. However, diﬀerent sets
of parameters very well could complicate the error surface providing a better use case for
parallel tempering.

While the need to balance the exploration-exploitation trade oﬀ was not prevalent in
this work, we still prefer the multi-replica algorithm over the single replica one because our
algorithm scales with negligible computational cost as replicas are added. It is worth noting
that care should be taken in the number and spacing of replicas [190]. Additionally, our
choice for the optimal β was not the best performing replica as it intentionally allows for
missteps to avoid local optima. However, without running multiple long simulations it is
not possible to know the complexity of the error surface. Therefore, we suggest using our
method for determining an optimal β and ensuring replicas with higher and lower values are
added to the simulation.

This work shows that it is possible to miniaturize a network, while preserving multiple
attributes. Future studies should investigate how the reduction in the graph relates to the
error. For example, if the graph is miniaturized by too large of an amount, it will no

71

longer be possible to match metrics, and conversely a small reduction should prove easy
to match metrics. However, the exact nature of this relation is unknown. We examined
a set of metrics we thought were important, but there is a quasi-inﬁnite number of metric
combinations that could be studied. At a minimum, exploring various combinations of
metrics to determine the diﬃcult in matching them would be useful for future studies. In
a similar vein, the initialization of graphs can be done such that they start closer to the
target metrics. Perhaps the use of greedy algorithms before the parallel tempering could
aid in ﬁnding optimal solutions. Most importantly, future work should study the amount of
error incurred from the miniaturization process and tie these errors to the metrics that are
matched. Speciﬁcally, investigating how much a result should be expected to change as a
function of reduction amount and metrics preserved.

The work presented here is an important ﬁrst step to miniaturizing social network for
use in simulations. Our parallel tempering algorithm with its excellent scaling and ability
to match metrics is a valuable contribution to optimizing graphs.

72

CHAPTER 6:

STOCHASTIC DIFFERENTIAL EQUATION BASED

MODELING OF DETERRENCE

6.1 Summary
This chapter begins the shift to the ﬁnal facet: machine learning. The main objective through
this chapter the development of a diﬀerential equation-based model of deterrence, the results
of which will inform decisions in future work on applying ABMs to modeling conﬂict. The
development of this model is my main contribution in the chapter. The analysis of this model
is in its early stages. Nonetheless, two types of global sensitivity analyses are performed on
a version of the model geared towards modeling short-term conﬂicts. Utilizing the Python
package SALib a variance-based Sobol and a PAWN global sensitivity analysis are performed
to examine the sensitivity of conﬂict outcomes to various parameter of the model. Two
proﬁles for competing actors are deﬁned and the sensitivity analyses are performed on the 3
possible combinations of actors.

Introduction

6.2
Conﬂict and deterrence have been central themes in the study of international relations and
military strategy for centuries [191, 192, 193]. The relevance of these themes has evolved
with the changing nature of warfare, prompting the need for sophisticated analytical tools
to understand and predict conﬂict dynamics. The use of mathematical and computational
methods in this domain can be traced back to pioneers like Thomas Schelling, whose work
in game theory and strategic analysis laid the groundwork for quantitative approaches to
conﬂict studies [194]. The importance of such studies lies not only in their academic value but
also in their practical applications. Governments and military strategists have long sought
eﬀective means to deter potential adversaries, and understanding the underlying mechanics
of conﬂict can inform more eﬀective policies and strategies. The mathematical modeling of
conﬂict allows for the analysis of complex systems and interactions, providing insights that
are often not apparent through qualitative methods alone.

Over the years, various models have been developed to analyze conﬂict, each with its
unique focus and limitations. For instance, the Lanchester model [195], is renowned for its
application in conventional warfare scenarios, providing insights into force attrition and the
outcomes of direct engagements. Many real world battle have been studied through the lens
of Lanchester [196, 197, 198, 199], and various extensions of Lanchester models have been
proposed, such three-way combat [200], irregular warfare [201], and adding spatial compo-
nents [202]. Another example, the Richardson model [203], delves into the dynamics of arms

73

races, highlighting the feedback loops inherent in competitive armament policies. This model
has also applied to Vietnam war [204] and middle east [205]. Guerrilla warfare models, such
as the one developed by Deitchman [206], oﬀer perspectives on asymmetric conﬂicts where
conventional military metrics may not apply. While these models have signiﬁcantly con-
tributed to our understanding of conﬂict dynamics, they primarily address tactical aspects
of warfare. Strategic, long-term considerations, particularly in the realm of deterrence, have
been less emphasized. Additionally, these models often do not account for the psycholog-
ical elements of conﬂict, such as the mindset of actors regarding deterrence, nor do they
adequately incorporate the stochastic nature of long-term conﬂict dynamics.

This chapter introduces a mathematical model that extends the traditional framework
of conﬂict analysis to address these gaps. Our model integrates the foundational principles
of the Lanchester model with new components that reﬂect the strategic, long-term aspects
of conﬂict and deterrence. Key innovations of our model include:

1. A system of ordinary diﬀerential equations (ODEs) that goes beyond ﬁght-to-the-ﬁnish

scenarios, reﬂecting long-term strategic considerations.

2. Explicit accounting for the role of deterrence, not just as a consequence of initial

conditions but as an integral component of the conﬂict dynamics.

3. Inclusion of psychological state equations, modeled as stochastic diﬀerential equations
(SDEs), to represent the deterred or undeterred mindset of actors, incorporating a
pitchfork bifurcation mechanism.

4. The introduction of noise terms in the psychological state equations to capture random
ﬂuctuations and events over extended periods. Coupling terms that link the psycho-
logical state of actors to their resource levels, aligning with the deterrence model.

5. Additional terms to enhance realism, including a surrender term that allows actors
to withdraw under certain conditions and a refresh term that simulates rebuilding
capabilities between conﬂicts.

By incorporating these elements, our model provides a more comprehensive and realistic
framework for analyzing conﬂict and deterrence over extended periods. This approach not
only ﬁlls existing gaps in the literature but also oﬀers practical insights for policymakers and
strategists seeking to understand and navigate the complex landscape of modern conﬂicts.

74

6.3 Methods
Our model is detailed in this section. Because the model describes several processes, we
develop the model progressively in several subsections beginning with a review of the Lanch-
ester model in the next subsection. We extend this model to include a resource refresh term.
In the following subsection we introduce a psychology model that reﬂects the mindset of the
actors; this model is stochastic to reﬂect long-term random events that trigger conﬂicts and
couples to the extended Lanchester model.

6.3.1 Lanchester’s models of conﬂict: Review
Lanchester models are a widely used class of conﬂict models used to study combat dynamics.
Typically consisting of systems of coupled ordinary diﬀerential equations, these models pre-
dict the outcome of ﬁght-to-the-ﬁnish conﬂicts. The models revolve around key parameters,
αx and αy, representing the combat strengths of two opposing actors with resource levels x
and y. For simplicity, we label the actors by the resources they control, i.e. actor x and
actor y. While Lanchester models often look similar, they model diﬀerent quantities, e.g.,
personnel or armaments. Here, we treat x and y as general resources that can be adapted
to other quantities through post-processing. This ﬂexibility allows for broader applications,
such as equating resource depletion to speciﬁc outcomes like causalities in a battle. Minguel-
Castro et al.
[193] overview various adaptations of Lanchester models tailored to diﬀerent
conﬂict scenarios. However, we build upon Lanchester’s foundational models [195], aiming
to integrate stochastic events and psychological aspects of deterrence, thereby extending the
models’ applicability to modern, multifaceted conﬂict environments.

Lanchester proposed two fundamental models to predict conﬂict dynamics in one-on-one
and all-on-all scenarios. Both of these models describe the time evolution of the number of
combatants available to each actor. The one-on-one model, formulated as

dx
dt
dy
dt

= −αyyx,

= −αxxy,

(6.1)

(6.2)

assumes the loss rate of each actor is proportional to the size of the opposing force and their
own force. Such a model resembles scenarios akin to series of individual duels among the
combatants. In contrast, the all-on-all model represented by

dx
dt
dy
dt

= −αyy,

= −αxx,

75

(6.3)

(6.4)

assumes losses are proportional to the size of the opposing force resembling a scenario where
combatants are engaged in collective battle. This fundamental distinction in the models
provides diﬀerent insights into conﬂict dynamics, with the one-on-one model emphasizing
the impact of mutual engagement, while the all-on-all model focuses on the sheer force size.

Figure 6.1: Example of Lanchester’s laws. In this example each actor’s available resources are
displayed with either a solid or dashed lines. The two colors represent which of Lancheser’s
laws describe the dynamics of the actors’ resources.

A comparison of these models is illustrated in Fig. 6.1, with parameters set at αx = 4
and αy = 8. The initial resources available to actor x are depicted in blue, and those of actor
y in orange. Solid lines represent the one-on-one model, while dashed lines correspond to
the all-on-all model. In both scenarios, actor x emerges as the victor. Notably, the all-on-all
enters a nonphysical regime once an actor’s resources are depleted; traditionally, the results
past this point are disregarded. Conversely, the one-on-one is unaﬀected by this issue due
to its inherent non-linear dynamics. This example, among many possible parameter sets
and initial conditions, illustrates the foundational principles of these models. However, from
these models general laws can be derived. From Eq. 6.1 and Eq. 6.2 Lanchester’s linear
law can be derived. This law states that team x wins if αxx0 > αyy0. That is, there is an
equal exchange for the number of ﬁghters and their strength in determining the outcome.
Similarly from Eq. 6.3 and Eq. 6.4 Lanchester’s square law can be derived. This law states
that team x wins if αxx2
0. Here the initial number of ﬁghters have a disproportional
importance compared to their strength.

0 > αyy2

Due to the non-physical behavior of the all-on-all model, we will use only the one-on-one

76

variant of Eqns. 6.3 and 6.4 in what follows. Additionally, when ﬁtting Lanchester’s models
to real battle data, researcher have found similar performance among the all-on-all and one-
one-one models [197]. In the next subsection we add a term needed to extend the Lanchester
model to longer time scales over which actors can regroup and replenish their resources. We
then turn to adding the psychological state of the actors x and y.

6.3.2 Lanchester’s models of conﬂict: Replenish
We are interested in long-term strategic modeling, and thus we require a term that allows
actors to replenish their resources. We introduce a resource equation of the form

dx/dt = g(1 − x/C),

(6.5)

which includes a growth rate g of resources and a term (1 − x/C) that limits the growth to
a capacity C. With this term alone, the resources obey

x = C − (C − x(t0))eg(t0−t)/C,

(6.6)

which steadily tends toward C on time scale g−1. Combined with the Lanchester model, we
obtain

dx
dt
dy
dt

= −αyyx + gx (1 − x/Cx) ,

= −αxxy + gy (1 − y/Cy) .

(6.7)

(6.8)

6.3.3 Psychology model
One of the main modiﬁcations of Lanchester models we considered was the addition of a
psychological state for the actors. We incorporated such a state as a random walk on a
potential energy surface. The stochastic diﬀerential equations that dictate the psychological
state consist of multiple terms that capture physical phenomenon. In Section 6.3.3 we discuss
the ﬁrst of such terms that constructs a quasi-binary state resembling the willingness of an
actor to engage in conﬂict. Then in Section 6.3.3 we discuss the threat-response term that
tries to match the psychological states of the actors. Next in Section 6.3.3 we discuss a term
that captures an actor surrendering due to excessive losses. Lastly, we address our deterrence
in Section 6.3.3 that utilizes Lanchester’s square law to inform an actor’s psychological state.

77

Quasi-binary state

To introduce a quasi-binary state into our model, we utilized a pitchfork bifurcation, com-
monly characterized by two stable and one unstable ﬁxed points under certain conditions.
We deﬁne the psychological potential function as

U (S) = S4 − 2S2,

(6.9)

where S represents the psychological state of an actor. This potential, visualized in Fig.
6.2, has minima at S = ±1 and a local maximum at S = 0. These minima correspond to
distinct psychological states: near S = −1 is the soft war state, characterizing a period of
non-engagement, and near S = 1 is the total war state, representing full combat engagement.
The stability of these minima and instability of the local maximum at S = 0 imply that
an actor’s state will settle near one of these two extremes. The variable γ controls the
randomness of transitions between the soft and total war states. Larger values of γ allow
for more frequent shifts between soft and total war, reﬂecting the unpredictable nature of
an actor’s decision-making process under varying external pressures

Figure 6.2: Psychological state potential energy surface. When an agent is in the left well,
they are in a peacetime state, and utilize few resources. If an actor is in the right state, they
will use resources to their fullest extent. Using this potential to describe the psychological
state of the actors creates a quasi-binary state system.

78

An actor’s psychological state is described by the stochastic diﬀerential equation

dS
dt

= −4S3 + 4S + γη,

(6.10)

where η represents non-correlated Gaussian white noise, characterized by (cid:104)η(cid:105) = 0 and
(cid:104)η(t)η(t(cid:48))(cid:105) = δ(t − t(cid:48)). This equation captures the dynamics of the psychological state
(S), where the term −4S3 + 4S reﬂects the tendency of the state towards soft and total war.
The addition of γη simulates the unpredictable nature of psychological responses in conﬂict
situations.

The evolution of this psychological state is illuminated through the Fokker-Planck equa-

tion, which describes the evolution of the probability density function P (S, t):

∂P
∂t

=

∂
∂S

(cid:18) dU
dS

(cid:19)

P

+ γ

∂2P
∂S2 ;

(6.11)

(see [207] for a review). This equation predicts how the distribution of psychological states
evolves, providing insights into the likelihood of diﬀerent states over time. The steady-state
probability density function is given by

P (S) =

N
γ

(cid:18)

exp

−

(cid:19)

,

U (S)
γ

(6.12)

where N is a normalization constant. This equation describes a long-term view of the
psychological states, indicating the probabilities of an actor being in various states under
stable conditions.

In Fig. 6.3 we visualize the psychological potential (shown in blue) with the steady-state
probability density function P (S) for diﬀerent values of γ – .25 in orange, 1 in green, and 4
in red. This illustration highlights how the size of random steps that is controlled through γ
inﬂuences the likelihood of an actor’s psychological state transitioning between states. When
the size of random steps is lower, indicated by lower values of γ, the psychological state is
most likely to be near one of the minima, as indicated by the pronounced peaks in P (S).
Conversely, when γ is higher, the size of random steps is larger and an actor’s psychological
state is less inﬂuenced by the potential. This is indicated by P (S) tending towards a uniform
distribution. Looking ahead, we will further develop our model by assigning individual
psychological states Sx and Sy to each actor, laying the groundwork for additional, more
complex terms.

79

Figure 6.3: Solutions of Eq. 6.12 for various values of γ. The psychological potential is
displayed in blue, with three values of the size of random steps in Eq. 6.10. The orange,
green, and red lines display the probability of a random walker being found at a certain
location. When the step sizes are small, the random walker is most likely to be found
in one of the wells. As the temperature increases, the distribution approaches a uniform
distribution, meaning the random walker is likely to be found in or between the wells.

Threat-Response

Building upon our model’s framework of psychological states, we now introduce a new com-
ponent: the threat-response term. This term is designed to dynamically adjust an actor’s
psychological state in response to perceived changes in their enemy’s stance. When an actor
perceives that their enemy is in a higher psychological state, this term acts to elevate the
actor’s own state, potentially increasing the rate resources are utilized. Conversely, if an
actor’s psychological state is higher than their enemy’s, the term works to moderate their
state, aligning closer to the enemy’s level. This mechanism allows for modelling strategic
adjustments made by actors in response to their enemy’s actions. The relative strength
of the threat-response term is controlled by a constant β; higher values of β imply more
pronounced adjustments.

For actor x, the threat-response term is expressed as β(Sy − Sx), and for actor y it is

80

β(Sx − Sy). We integrate the threat-response term with the quasi-binary state as

dSx
dt
dSy
dt

= −4S3

x + 4Sx +

= −4S3

y + 4Sy +

β
2
β
2

(Sy − Sx),

(Sx − Sy),

(6.13)

(6.14)

where each actor has their own psychological state denoted by a subscript. To gain insight
on how this term alters the dynamics, we examine the location and stability of ﬁxed points
of the system, listed in Table 6.1 and visualized in Fig. 6.4.

Sx
0
1
−1
√
4 − β/2
√
−
4 − β/2
√
−Z − β + 8/4
√
−Z − β + 8/4
√
Z − β + 8/4
√
Z − β + 8/4

−

−

Sy
0
1
−1

(4 − β)(
−(4 − β)(

√
4 − β − 4 + β/2)/β
√
4 − β − 4 + β/2)/β

√
8 − Z − β(β − Z − 8)/8
√
8 − Z − β(β − Z − 8)/8
−
√
Z − β + 8(β + Z − 8)/8
√
Z − β + 8(β + Z − 8)/8

−

Table 6.1: Fixed points and stability of ﬁxed points for threat-response term. Here Z =
(cid:112)−3β2 − 16β + 64.

In our analysis, nine ﬁxed points are identiﬁed, though their presence and positions
depend on the values of β. Particularly noteworthy are β = 0, β = 8/3, and β = 4, where
the ﬁxed points tend to align or overlap.

These these alignments and convergences, are demonstrated in Fig. 6.4 through quiver
plots for β = 0.01 (Fig. 6.4a), β = 1 (Fig. 6.4b), β = 3 (Fig. 6.4c), and β = 5 (Fig. 6.4d). In
these ﬁgures, stable nodes, saddle points, and unstable nodes are represented as blue squares,
purple circles, are red diamonds, respectively. The arrows indicate the derivatives of Sx and
Sy with respect to time, illustrating the evolution of the actor’s psychological states and the
speed of change.

When β is near zero, the system has nine ﬁxed points: four stable nodes at (±1, ±1) and
(±1, ∓1), four saddles points at (0, ±1) and (±1, 0), and an unstable node at (0, 0). In this
regime, each actor acts independently, inﬂuenced predominately by the quasi-binary state.
This independence results in actors gravitating towards the nearest minima unless starting at
a psychological state of zero. As β increases, the saddle-points at (−1, 0) and (0, 1) approach
the stable node at (−1, 1), while the saddle points (0, −1) and (1, 0) approach the stable
node at (1, −1). Notably these groups of ﬁxed points are approaching (0, 0) along the line

81

Sx = −Sy. At β = 8/3, these points collide, transforming into two saddle points. With
further increases to β, these saddle points continue to move towards the central unstable
node, until at β = 4, they merge into a single saddle point at (0, 0). Throughout these
changes, the stable nodes at (±1, ±1) remain ﬁxed. When β exceeds 4, actors are more
inclined to align their psychological states with their enemies. If positioned symmetrically
on either side of zero, (on the line Sx = −Sy), actors are drawn to the (0, 0) state. Otherwise,
their states are attracted to a soft war for negative dominate psychological states, or total
war for positive ones.

The dynamics of these ﬁxed points align with the anticipated behavior of a defensive
term in our model. Actors who assign a lower weight to their enemy’s psychological state
remain largely unaﬀected by their enemy’s state, indicating independent decision making.
In contrast, as the weight increases, actors tend to mirror their enemy’s psychological state,
reﬂecting adaptive strategies in response to changing conﬂict dynamics. Crucially, all ﬁxed
points are bounded within the region [−1, −1] × [1, 1], ensuring that the values of Sx and
Sy remain in this realistic range, regardless of the value of β. This boundedness prevents
the psychological state of the actors from reaching unrealistic extremes such as ∞ or −∞,
thereby maintaining the model’s applicability to real-world scenarios.

Surrender

Having established the dynamics of the threat-response term, we now turn our attention
to another aspect of our model: the surrender term. This term plays a role in simulating
the conditions and consequences of surrender, and further enhancing our understanding of
decision-making in conﬂict scenarios. Comprising of two components, the surrender term
ﬁrstly acts as a switch that turns of the term once an actor’s psychological state reaches
S = −1. The second component models the reduction in an actor’s psychological state pro-
portional to a losses in resources, characterized by R. This captures the diminishing morale
and willingness to continue conﬂict as losses mount. The inﬂuence of this term is controlled
(cid:1), for
by the constant, (cid:15). For actor x, the term is (cid:15) (Sx + 1) (cid:0)x − 1
actor y.

(cid:1), and (cid:15) (Sy + 1) (cid:0)y − 1

R

R

Considering the coupling of these terms to the available resources, identifying and vi-
sualizing ﬁxed points for this system becomes challenging. Nonetheless, we can examine
the inﬂuence of the surrender term on the psychological state of an individual actor. For
simplicity, let’s consider the psychological state of actor x describe by the equation:

dSx
dt

= −(cid:15) (Sx + 1) (x − R) .

(6.15)

82

(a) When β approaches 0, the saddle-points align
with the stable points. This alignment in combi-
nation with the unstable node push the psycho-
logical state towards the nearest corner.

(b) Here we let β = 1. Once β hits 8/3 these
saddle-points combine with the stable nodes, re-
sulting in two saddle-points.

(c) Here we let β = 3. Once β is greater than
8/3, the resulting saddle-points tend towards the
center.

(d) Here we let β = 5. Once β passes 4 there is
a saddle-point at the center with stable nodes at
(±1,±1). Actors will either be in a total war or
peace state.

Figure 6.4: Stable points as a function of β, where blue squares are stable nodes, purple
circles are saddle points, and red diamonds are unstable nodes.

83

We simplify our analysis by assuming that losses, l = (x − R), remain constant. Using this
assumption, the solution to Eq. 6.15 is

s(t) = (s0 + 1)e−(cid:15)lt − 1,

(6.16)

where s0 is the initial psychological state of the actor. This result shows that the psycho-
logical state of the actor exponentially decays to −1, indicating an increasing inclination to
surrender as losses continue. The rate of this decay is controlled by the weight (cid:15) and losses, l.
An actor with a lower value of (cid:15) will have a slower reaction to accumulating losses, reﬂecting
a more resilient or less responsive attitude towards surrender, compared to an actor with a
higher (cid:15), who will react more swiftly. This variation can be interpreted as diﬀerent tolerances
for surrendering among actors in conﬂict scenarios.

Conﬁdence-Deterrence term

The ﬁnal component we introduce to our model is the integration of conﬁdence and deter-
rence factors. Traditional one-shot ﬁght-to-the-ﬁnish models implicitly incorporate deter-
rence within their logical structure, but often do not explicitly model it. An actor’s decision
to engage in battle is based on their perceived performance, e.g. evaluating one of Lanch-
ester’s laws. Only if an actor decides to engage is the model evaluated. In contrast, our
model seeks to explicitly incorporate deterrence into the psychological states of the actors.
This includes not only a weaker actor being deterred, but also the strong actor gaining
conﬁdence based on strategic assessments. Speciﬁcally, each actor evaluates their likelihood
of success by applying Lanchester’s square law, comparing the combat strengths αyy2 and
αxx2 weighted by δ. This approach allows us to simulate not only the physical aspects of
conﬂict but also the psychological warfare that often precedes and dictates the course of
engagements.

In our model, each actor evaluates their combat strength by subtracting it from their
enemy’s strength. For actor x this has the form αxx2 − αyy2 and for actor y its αyy2 − αxx2.
An actor who perceives a higher combat strength gains conﬁdence, leading to an increase
in their psychological state. Conversely, an actor who determines that they are at a combat
disadvantage is deterred, and their psychological state is decreased correspondingly.

6.3.4 Lanchester with deterrence
Using the discussed terms, we move to couple the equations for resources available to and
psychological states of the actors. The model is further generalized to allow for each actor

84

to separately weight each term, denoted by subscripts. The full model is given by:

dx
dt
dy
dt
dSx
dt
dSy
dt

= −αyyxB(Sy) + gx

1 −

(cid:18)

= −αxxyB(Sx) + gy

1 −

(cid:18)

(cid:19)

(cid:19)

,

,

x
Cx
y
Cy

= −4S3

x + 4Sx +

= −4S3

y + 4Sy +

βx
2
βy
2

(Sy − Sx) +

(Sx − Sy) +

(cid:15)x
2
(cid:15)y
2

(6.17)

(6.18)

(Sx + 1) (x − Rx) + δx(αxx2 − αyy2) + γxηx, (6.19)

(Sy + 1) (y − Ry) + δy(αyy2 − αxx2) + γyηy,

(6.20)

where B bounds Sy and Sx between −1 and 1. Each of the parameters and the range of
values they can take is listed in Table 6.2.

variable
0 ≤ x ≤ 1
0 ≤ y ≤ 1
−∞ ≤ Sx ≤ ∞
−∞ ≤ Sy ≤ ∞
0 ≤ βx(y) ≤ 1
0 ≤ (cid:15)x(y) ≤ 1
0 ≤ Rx(y) ≤ 1
0 ≤ δx(y) ≤ 1
0 ≤ γx(y) ≤ 1
0 < Cx(y) ≤ 1
0 < αx(y) ≤ 1
0 ≤ gx(y) ≤ 1

description
fraction of total available resources for actor x
fraction of total available resources for actor y
psychological state of actor x
psychological state of actor y
strength of threat-response term
strength of surrender term
losses an actor is willing to take willing to take
strength of conﬁdence-deterrence term
controls size of random jumps
maximum number of resources for an actor
ﬁghting strength of an actor
growth rate of an actor’s resources

Table 6.2: Model parameters and acceptable values.

While these parameters can be set to any number in their range, we construct a few sets
of parameters to represent types of actors. The sets are broken down into two types based
on an actor’s risk aversion. A conservative actor is less likely to accept losses, more likely to
surrender, have medium-level of conﬁdence, and more likely to respond to an enemy with a
matched response compared to a reckless actors. These traits are captured for conservative
actors by assuming R = 1, (cid:15) = 1, δ = .5, and β = 1, while for reckless actors we assume
In both cases, we set C = 1 to represent each actor
R = 0, (cid:15) = 0, δ = 1, and β = 0.
attempting to build resources to the same highest level.

Using our model and actor proﬁles, we performed sensitivity analyses on our model
for confrontations between diﬀerent actor proﬁles. We solved out model using the RK-45
implementation in scipy [208], and examined the sensitivity of the steady-state value of x, y,
Sx and Sy to the 4 remaining actor speciﬁc parameters and the initial conditions of resources

85

and psychological states. In practice we solved the system of equations until a simulated
time t = 100, with a max time step of .1. Utilizing the Python package SALib [209, 210], we
generated 1,024 Sobol quasi-random sequences, each sequence consists of 10 sets of values for
the model’s parameters and initial condition. In total, our model was executed 10,240 times.
Using these sets of inputs and corresponding outputs, we performed a Sobol, variance-based
global sensitivity analysis [104], and the moment-independent global sensitivity analyses
PAWN [211].

A variance-based global sensitivity analysis relates the variance in a model’s parameters
to variance its outputs. Typically this relation is presented with the ﬁrst-order (S1) and
total eﬀect (ST ) Sobol indices. The ﬁrst-order eﬀect measures how much of the model’s
output variance can be described by a single input on its own, while the total eﬀect captures
interactions among parameters as well. While often informative, variance-based methods are
not well suited for highly-skewed or bimodal output distributions. In these cases moment-
independent methods better capture the importance of inputs on outputs. In this study, we
utilized the PAWN method that measures the change in the cumulative distribution function
of the output when all inputs are allows to vary versus when inputs are held constant. This
change is measured using the Kolmogorov–Smirnov statistic.

6.4 Results
In this section we discuss the results of our sensitivity analyses conﬂicts between two con-
servative actors (Table 6.3), two reckless actors (Table 6.4), and a conservative and reckless
actor (Table 6.5). These results are presented in three tables, one for each pair of actors.
Each table is divided into four subtables where the ﬁnal states of the resources, denoted
In
x(∞) and y(∞), and psychological states, denoted Sx(∞) and Sy(∞), are presented.
each subtable, we bold the sensitivity index that corresponds to the parameter that each
model output is most sensitive.

In all scenarios, the distributions of values for x(∞) and y(∞) are left skewed with a peak
near 1, while the distributions of values for Sx(∞) and Sy(∞) are both bimodal, with peaks
at ±1 for the 10,240 simulations we examined. We performed a variance-based sensitivity
analysis on these simulations, but since the distributions of the ﬁnal states are not well
described by variance alone, we also present a moment-free sensitivity analysis.

In all of our analyses we see that the variance-based and moment-free analyses produce
similar rankings for which parameters x(∞), y(∞), Sx(∞) and Sy(∞) are sensitive. While
the ranking of parameters are similar, the diﬀerence in sensitivity indices is smaller in the
moment-free analysis. As expected, this is more evident in the sensitivities of the psycho-
logical states that have bimodal distributions, and thus are not well described by variance

86

alone. Nonetheless, all model outputs are most sensitive to the initial conditions of the
psychological states, Sx(0) and Sy(0), with x(∞) and Sy(∞) most sensitive to Sy(0), while
y(∞) and Sx(∞) most sensitive to Sx(0). The observation that the ﬁnal values of the psy-
chological state are most sensitive to their initial conditions agrees with the behavior we
see in Fig. 6.4. This ﬁgure examines the behavior of the psychological state when only the
quasi-binary state and threat-response term are at play. Speciﬁcally, in Fig. 6.4a when the
threat-response term is small there are four stable nodes that determine the possible ﬁnal
values of the psychological states that are completely controlled by the initial conditions of
the psychological states. Even as the threat-response term becomes larger there are consis-
tently at least two stable nodes that determine the ﬁnal values of the psychological state,
and which value the model tends towards is fully determined by the initial value of the psy-
chological states. While there are other terms at play, they must not be strong enough to
overcome the eﬀects of the quasi-binary state term.

A reckless actor’s ﬁnal psychological state has a low sensitivity to their opponent’s initial
psychological state this is partly expected as when (cid:15) = 0, some inﬂuence of the opponent’s
psychological state is removed, but nonlinear eﬀects can still enter through the conﬁdence-
deterrence term. However, a conservative actor is inﬂuenced by their opponent’s initial
psychological state. When examining the PAWN indices, conservative actors are about twice
as sensitive to their own initial psychological state as their opponent’s initial psychological
state, while reckless actors are over 6 times more sensitive to their own initial psychological
state compared to their opponents.

The ﬁnal value of an actor’s resources are second most sensitive to growth rate of their
resources and their enemy’s ﬁghting strength. A reckless actor is about slightly more sensitive
to their own ﬁghting strength compared to a conservative actor. Additionally reckless actors
are more sensitive to their own initial conditions compared to conservative actors.

87

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

T
S1
ST
0.030
0.000
0.010
0.036
0.002
0.043
0.105
0.160
0.015
0.658 0.387 0.361
0.116
0.013
0.033
0.138
0.109
0.213
0.139
0.169
0.275
0.053
0.005
0.049

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

(a) x(∞)
S1
ST
0.010
0.074
0.019 −0.003
0.803
0.976
0.016
0.142
0.004
0.038
0.046
0.000
0.026 −0.001
0.003
0.007

(c) Sx(∞)

T
0.039
0.028
0.454
0.268
0.117
0.113
0.064
0.049

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

S1
ST
0.011
0.061
0.024 −0.005
0.676
0.404
0.113 −0.022
0.105
0.195
0.007
0.054
0.004
0.052
0.180
0.265

(b) y(∞)
S1
ST
0.000
0.011
0.012
0.061
0.038
0.147
0.820
0.962
0.009
0.034
0.000
0.022
0.002
0.002
0.016 −0.007

(d) Sy(∞)

T
0.039
0.031
0.357
0.156
0.151
0.115
0.050
0.133

T
0.027
0.034
0.264
0.453
0.121
0.104
0.042
0.068

Table 6.3: Global sensitivity analysis for conservative actor vs conservative actor. Here ST
is the total eﬀect Sobol index, S1 is the ﬁrst-order eﬀect Sobol index, and T is the PAWN
index. For convenience, we bold the most important input as judged by each index. The
ﬁnal value of resources is most sensitive to their opponent’s initial psychological state, while
the the ﬁnal value of the psychological state is most sensitive to its own initial condition.

88

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

S1
ST
0.010
0.037
0.046 −0.009
0.028
0.053
0.437
0.643
0.022
0.057
0.137
0.247
0.164
0.252
0.003
0.029

(a) x(∞)
S1
ST
0.005
0.043
0.000
0.043
0.894
0.988
0.000
0.000
0.006
0.056
0.020
0.058
0.012
0.000
0.006 −0.002

(c) Sx(∞)

T
0.028
0.029
0.074
0.413
0.132
0.148
0.208
0.042

T
0.030
0.027
0.500
0.078
0.197
0.200
0.041
0.058

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

T
S1
ST
0.032
0.008
0.041
0.033
0.035
0.002
0.675 0.453 0.407
0.074
0.002
0.041
0.168
0.136
0.236
0.140
0.027
0.070
0.041
0.004
0.032
0.216
0.150
0.250

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

(b) y(∞)
S1
ST
0.016
0.035
0.054 −0.011
0.000
0.000
0.892
0.990
0.011
0.056
0.042
0.004
0.002 −0.002
0.000
0.000

(d) Sy(∞)

T
0.029
0.031
0.076
0.500
0.203
0.201
0.053
0.049

Table 6.4: Global sensitivity analysis for reckless actor vs reckless actor. Here ST is the total
eﬀect Sobol index, S1 is the ﬁrst-order eﬀect Sobol index, and T is the PAWN index. For
convenience, we bold the most important input as judged by each index. The ﬁnal value of
resources is most sensitive to their opponent’s initial psychological state, while the the ﬁnal
value of the psychological state is most sensitive to its own initial condition.

89

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

T
S1
ST
0.034
0.002
0.012
0.044
0.002
0.051
0.121
0.170
0.023
0.638 0.378 0.354
0.109
0.010
0.032
0.142
0.111
0.218
0.140
0.180
0.284
0.054
0.009
0.050

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

(a) x(∞)
S1
ST
0.005
0.043
0.000
0.042
0.894
0.987
0.000
0.000
0.007
0.056
0.020
0.059
0.012
0.000
0.006 −0.002

(c) Sx(∞)

T
0.031
0.029
0.500
0.073
0.190
0.186
0.043
0.059

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

parameter
x(0)
y(0)
Sx(0)
Sy(0)
αx
αy
gx
gy

S1
ST
0.007
0.040
0.000
0.032
0.651
0.429
0.028 −0.003
0.135
0.235
0.031
0.071
0.005
0.040
0.165
0.265

(b) y(∞)
S1
ST
0.000
0.011
0.007
0.073
0.047
0.170
0.811
0.962
0.027 −0.005
0.030 −0.007
0.007
0.000
0.029 −0.005

(d) Sy(∞)

T
0.033
0.036
0.403
0.066
0.167
0.142
0.046
0.208

T
0.026
0.034
0.304
0.455
0.121
0.113
0.039
0.067

Table 6.5: Global sensitivity analysis for reckless actor (x) vs conservative actor (y). Here
ST is the total eﬀect Sobol index, S1 is the ﬁrst-order eﬀect Sobol index, and T is the PAWN
index. For convenience, we bold the most important input as judged by each index. The
ﬁnal value of resources is most sensitive to their opponent’s initial psychological state, while
the the ﬁnal value of the psychological state is most sensitive to its own initial condition.

90

6.5 Conclusions
In this work we proposed a model for long-term strategic studies of the statistical properties
of conﬂict. This model was adapted from a Lanchester model with replenishment that was
further coupled to a psychology model. Our psychology model included multiple terms
representing responding to threats, surrendering, conﬁdence, and deterrence. We deﬁned
multiple proﬁles for actors involved in a conﬂict, and performed extensive sensitivity analyses
to examine the outcomes of matches between these actor proﬁles. These results informed the
importance of ﬁghting capabilities, resource levels, and psychological states on the outcome
of the conﬂict.

We examined two types of sensitivity analyses based on the distributions of outputs from
our model that produced similar rankings among the parameters. Overall, we saw that the
initial conditions of the psychological state were the dominate parameter in determining
the ﬁnal state of the model in all scenarios. This result was not unexpected based on the
examinations of subsets of the model terms in Sec. 6.3. Even though the quasi-binary state
has values between 0 and 1, in future iterations of this model, an additional parameter should
control its relative strength to other terms. An alternative could be tuning the noise term
such that it allows for the psychological state to switch.

When examining the ﬁnal level of resources among the actors, we found this ﬁnal value
was second most sensitive to how quickly the resources are replenished and how strong the
opponent’s ﬁghting strength was. Such an observation implies that a stronger enemy could
be successfully defended by replenishing resources faster than they can destroy, even if your
ability to ﬁght is less than your opponents. For both types of actors this behavior was
present without a dependence on the proﬁle of the opponent, suggesting this observation is
robust to proﬁle type. When observing the ﬁrst-order eﬀect, we see much of the variance in
the amount of ﬁnal resources is being driven by interactions among the parameters. Future
work could examine higher order indicies that describe the eﬀects of subsets of parameters.
In our analysis, we also saw that reckless agent’s ﬁnal level of resources were more sensitive
to their initial conditions compared to conservative actors. This observation combined with
the one that reckless agent’s have a lower sensitivity to their opponent’s initial psychological
state than conservative actors suggest self-centered behaviors from the reckless agents. A
conservative agent who ﬁnds themselves in a conﬂict with such a reckless agent would beneﬁt
from adopting strategies that eﬀect an opponents resources before a conﬂict breaks out as
well as attempting to adjust their psychological states to more favorable values.

In future studies we hope to continue to explore theoretical conﬂicts using additional actor
proﬁles. These proﬁles can constrain the same parameters we have explored here, or even
be less restrictive. We note that in preliminary sensitivity analyses that allowed all model

91

parameters to vary resulted in result that were hard to map to actionable results. However,
performing such a large analysis did aid in model building by informing which parameters
the model was sensitive. Comparing this information to our expectations allowed for changes
to be implemented that were more aligned to our modeling goals. Additional actor proﬁles
can be mapped to real world leaders to help inform strategies in a conﬂict. The choices made
here were coarse, but a fair starting point.

While our analysis here focused on a deterministic version of the model, we plan to address
including stochastic terms in future work. The addition of stochastic terms will require new
quantities of interest to be examined. In preliminary runs of a stochastic model, we found
the number of, size of, and duration of conﬂicts were interesting and dynamic quantities.
Another topic worth exploring is the use of other base models. In this work, we extended
a two player version of Lanchester’s square law. However, there are many other versions,
such as adding a third player, that provide an opportunity for richer conﬂict dynamics. A
ﬁnal area of interest for future work is improvements to our psychological potential energy
surface. Here, we only included a static potential that captured a two state system. However,
adjusting the parameters of this surface, as well as making it time dependent could lead to
informative new insights. While the development of this model is in an early stage, it shows
promise for expanding the capabilities of ordinary diﬀerential equation based models to
explicitly account for deterrence. Additionally, insights gained from this model can lead to
improved experiments in more complex models, such as agent-based models, of conﬂict.

92

CHAPTER 7:

APPLYING REINFORCEMENT LEARNING TO

AGENT-BASED SIMULATIONS OF CONFLICT

7.1 Summary
In this chapter we shift fully to incorporating machine learning into agent-based models.
The main objective is to create an agent-based implementation of the game capture the
ﬂag such that reinforcement learning techniques can be applied to control players in the
game. The reinforcement learning packages Gymnasium [212] and stable-baselines3 [213]
are utilized to frame the implementation of capture the ﬂag as a reinforcement learning
problem and subsequently train a reinforcement learning agent. The reinforcement learning
agent is incrementally trained to solve increasingly complex tasks. These tasks starting with
controlling a single player to capture an undefended ﬂag and work towards controlling a
team to play against another team. Additionally, aspects of deterrence are examined in the
game. Namely, the willingness of a reinforcement learning agent to play as the penalty it
receives for being captured increases. While this work is preliminary, there is an extensive
discussion of paths to continue this project.

Introduction

7.2
Conﬂict is a complex phenomenon that is inﬂuenced by various factors ranging from indi-
vidual diﬀerences to societal contexts. For the purposes of this thesis, we deﬁne conﬂict as
a scenario where two or more groups have incompatible interests, and the actions of these
groups interfere with one another. On a personal level, conﬂict arises from diﬀerences in
relationships, values, facts, structures, and interests [214, 215]. While personal disputes
arise from these individual diﬀerences, on a macro scale, the dynamics shift and evolve. At
this broader scale, power dynamics, climate [216], inequality [217], and cultural and social
contexts play signiﬁcant roles [218]. Beyond these, broader material and existential reasons
further incite conﬂict [219]. Addressing and understanding conﬂict demands a holistic, mul-
tidisciplinary approach, merging insights from diverse ﬁelds including psychology, sociology,
political science, and computer science. This understanding can aid in anything from resolv-
ing personal disagreements to mediating national-scale disputes. A particularly intriguing
dimension of conﬂict is its asymmetric nature, where the entities in opposition have unequal
resources.

Asymmetric conﬂicts are studied in many contexts. For example, Clark and Konrad
explore multifront wars, illustrating scenarios where a defending team is spread across mul-
tiple fronts. In contrast, the attacking team is only required to capture a single front to be

93

victorious [220]. In another study, Dunne et al. propose a three-stage game where a weaker
challenger triumphs over a stronger incumbent, leveraging unconventional means that catch
the incumbent oﬀ-guard [221]. A distinct exploration of asymmetric conﬂict dynamics can be
seen in the works of Johnson and MacKay, who explore the implications of Lanchester’s laws
[222] in modern warfare [223]. These laws predict the outcome of conﬂicts, either one-on-one
or all-against-all, by evaluating the size and combat capabilities of the engaged groups. The
parameters of Lanchester’s laws make them particularly well-suited for studying asymmetric
conﬂict. Chapter 6 discusses extended Lanchester models with added replenishing terms
and a model of actor’s psychological states. Building on the extended Lanchester models,
this chapter employs the game of capture the ﬂag (CTF) as a representative model to delve
deeper into asymmetric conﬂict dynamics and deterrence.

Capture the ﬂag is a popular adversarial game that has been used to study conﬂict in
cyber security [224, 225, 226] and robotics [227, 228, 229, 230], among other ﬁelds. In the
physical world, the game is typically played in a large area, such as a ﬁeld or gymnasium,
and consists of two teams that compete to capture each other’s ﬂag, while simultaneously
defending their own. Gameplay unfolds through a series of strategic decisions, such as
whether to defend one’s own ﬂag or attack the opponent’s ﬂag, and how to coordinate
with one’s team members. Computational adaptations of CTF follow the general rules of
the game, however the inherent ﬂexibility of the medium allows for adaptations that may
be unfeasible to study in a physical game of the CTF. For example, in a computational
implementation players ability to detect their surroundings can be varied. Such discrepancies
can be mapped onto real conﬂict scenarios, and allow CTF to be used to study asymmetric
conﬂict.

Capture the ﬂag lends itself well to the modeling technique reinforcement learning (RL).
Reinforcement learning is often framed as a Markov decision process (MDP), and revolves
around two fundamental components: an agent, which serves as the learning entity, and
an environment, with which the agent interacts. The agent takes a sequence of actions
within the environment over a deﬁned number of iterations, T . At each iteration, t, the
state of the environment is represented by st ∈ S, where S represents the set of all possible
conﬁgurations. The agent observes the current state, st, and selects an action, at, from
the set of possible actions, A. Performing the action transitions the environment from st
to st+1, while the agent receives a reward rt+1 from the reward set, R. Notably, the state
may remain unchanged (st = st+1), and the reward can be positive, negative, or zero. The
goal of the agent is to learn a policy, π, that represents the series of action that maximize
the expected cumulative reward. Typically, rewards are discounted using a discount factor
γ ∈ [0, 1], to balance short-term gains against long-term goals. Thus, the total cumulative

94

reward at iteration t is computed as

Rt =

T −t−1
(cid:88)

i=0

γirt+i+1,

(7.1)

where an agent with γ close to 0 exhibits shortsighted behavior, while γ closer to 1 implies
a more farsighted approach.

Applied to CTF, a team is controlled by a reinforcement learning agent, similar to a
coach or general directing the actions of the team. The positions of players and ﬂags dictate
the state of the system. As dictated by the reinforcement agent, players take actions to
update their positions, and are correspondingly rewarded based on events in the game,
such as capturing the ﬂag or getting captured. There are many algorithms to solve RL
tasks [231, 232], many of which are performed in discrete action and state spaces. Such
as Q-learning [233], deep Q-learning [234], value-iteration [235], and temporal diﬀerence
[236]. Some of these methods can handle continuous state spaces, but fail when both the
state and action space are continuous. Other algorithms have been developed to solve such
problems [237, 238, 239, 240]. Recently genetic algorithms have been applied to training
reinforcement learning agents, and been found to be faster to train and produce comparable
solutions as other methods in certain scenarios [241, 242]. The Python package stable-
baselines3 [213] provides implementations of many standard RL algorithms. In particular,
their implementation of the Proximal Policy Optimization (PPO) algorithm is utilized in
this chapter.

Our overarching goal is to use reinforcement learning to train teams in our capture the ﬂag
environment that have various levels of ability. These abilities can include, how individuals
perceive and move through the environment, the size of teams, and the length of time a team
is trained. By competing these teams against each other, we can examine their deviations
from simpler models, such the Lanchester model discussed in Chapter 6. Such an exploration
can be used as a proxy to study asymmetric conﬂict in the real world.

7.3 Methods

7.3.1 Creating Capture the Flag Environment
There are many variations of CTF that vary in the number of players and ﬂags, as well as the
rules for capturing these entities. Traditionally, CTF is played with two teams, in a rectan-
gular, two-dimensional space such as a ﬁeld. The goal of each team is to capture the enemy’s
ﬂag and return in to their side without being captured while simultaneously defending their
own ﬂag. For our implementation, we have relaxed some of these assumptions.

We implemented a version of CTF in Python using Shapely [243] that allowed us to

95

easily deﬁne elements of the games as geometries and test for their intersections. We refer
to the area the game unfolds as the ﬁeld, and it is represented as an adjustable rectangle
Ω = [xmin, xmax] × [ymin, ymax] ∈ R2. The ﬁeld is split into a left and right side, with an
impenetrable boundary that keeps the players within its bounds. One team defends the
right side ((xmax − xmin)/2 > 0), while the other defends the left side ((xmax − xmin)/2 < 0).
Players and ﬂags are represented as circles in Ω with radii Rp and Rf , respectively. A team’s
players and their ﬂag initially start on the team’s side of the ﬁeld, but players can move
throughout the ﬁeld as the game progresses. The game is broken into a series of iterations,
where players select a velocity vector, v, to move themselves to a new location, with a
maximum magnitude of vmax. During each iteration there are two events that can occur,
depending on entities overlapping. These events are: 1) a player getting captured, when a
two opposing players overlap, 2) a ﬂag getting captured, when a player overlaps an opposing
ﬂag, and 3) a player hits the boundary, when a player overlaps with the boundary (see Table
7.1). When players overlap, the player that is on the side they are defending captures the
opposing player, removing them from the game. When a player overlaps with the opposing
team’s ﬂag, the ﬂag is captured and the game ends. If a player intersects the boundary,
they are moved back in bounds, perpendicular to the intersection. The team that is able
to capture the enemy’s ﬂag without having all their players captured wins. If neither team
successfully captures the enemies ﬂag within 1000 iterations, they both lose.

7.3.2 Applying Reinforcement Learning to Capture the Flag En-

vironment

With the basics of game play deﬁned, players require a method to determine how they
move during a game. One method would be to supply each player with a series of rules to
determine their next move, however we do not want to limit ourselves to strategies that we
propose. Instead, we wish to learn policies that can adapt to what other players are doing.
To accomplish this, we framed the game as a reinforcement learning task.

Each simulation consists of two team comprised of N1 and N2 players. The players on
each team are supplied actions either by an agent or a supplied set of rules. The players on
a team are unable to fully detect their environment, and instead rely on a set of R rays of
length L that are evenly spaced around their body. These rays can detect the distance from,
type of, and velocity of objects they intersect. The possible object types are teammate,
enemy, team ﬂag, enemy ﬂag, boundary, and open space. The agent (or rule set) that
controls each team can observe the positions and velocities of each team member, if the
member is defending or attacking, and the information detected by each ray. In total there
are (4R + 5)N1 and (4R + 5)N2 variables for the teams. As an example, a single iteration of a

96

simulation is shown in Fig. 7.1. Here, the bounds of the ﬁeld are given by xmin = ymin = −1
and xmax = ymax = 1, with the boundary represented as a black rectangle. The teams
are denoted with red and blue entities, where the larger circles are the teams ﬂags, and
the smaller circles are the individuals that make up the teams. The blue lines emanating
from the blue player are the rays it uses to perceive the environment. Most of the rays are
intersecting walls, however there are rays intersecting each ﬂag and the red player. Therefore,
the blue player can see these entities during this iteration. However, in a later iteration if
they rays are no longer intersecting these entities the blue player would not detect them. In
this example, only the blue players rays were shown to reduce clutter in the image.

Figure 7.1: Example capture the ﬂag environment. The boundary marking out-of-bounds
is shown in black. The larger colored circles are each team’s ﬂag, while the smaller circles
are the players. The rays emanating from the blue player allow the player to observe the
environment. Along each ray the player can detect what type of object the ray is intersecting,
the distance to the object, and the velocity the object is moving.

Using the input from the players on a team, an agent determines the action each player
should take. In this model, a player’s action is a velocity vector used to update its position.
The RL agent uses a neural network with two hidden layers, each with 140 neurons with
tanh activation functions to produce a velocity vector for each player. Each iteration after
each player has updated thier position, the RL agents receive rewards based on the player-
involved events that occurred. The events and their associated rewards are listed in Table 7.1.
Negative rewards are given in an eﬀort to penalize behaviors, while positive ones encourage
behaviors. In this model, each step taken incurs a reward of −1 to encourage the agent to

97

ﬁnish the game quickly. Similarly, a reward of −2 is given when players touch the boundary
to encourage them to stay in bounds. The ﬁnal negative rewards are given when a player
or their ﬂag is captured. Each of these rewards are relative to the number of iterations that
pass before they occur, such that they are 0 if the player or ﬂag are never captured in the
maximum number of iterations. Positive rewards are given for capturing the enemy ﬂag or
players. These rewards start at 1000 and are reduced by the number of iterations that pass
to encourage a quick solution.

event
player takes action
player touches boundary
player captures ﬂag
player captures enemy
player’s ﬂag is captured
player is captured

reward
−1
−2
1000 − iterations
1000 − iterations
iterations − 1000
iterations − 1000

Table 7.1: Rewards for capture the ﬂag. Negative rewards penalize behaviors while positive
ones encourage them. The ﬁrst two rewards encourage the RL agent to act quickly and keep
players in bounds. The second two rewards encourage the agent to capture entities, while
the ﬁnal two rewards encourage the agent to avoid being captured and protect their ﬂag.

Based on our implementation of CTF and associated rewards, we used the Proximal
Policy Optimization (PPO) algorithm implemented through stable-baselines3 [213] to train
RL agents to control a team. Our training occurred as a series of progressively more complex
tasks that are discussed below. In each of these simulations, we used xmin = ymin = −1 and
xmax = ymax = 1. The ﬂags had a radius of Rf = .1, while players had a radius of Rp = .05
and R = 16 rays of length of L = 2.

The simplest task we examined was training a single agent to capture a ﬂag without a
defender. For this scenario, the ﬁeld was not partitioned into sides, so the single player was
always in an attacking state. We made this choice to ensure that the agent examined the
ﬂag from all directions, and would not simply learn to move in a singular direction. Only
the ﬁrst three rewards in Table 7.1 are available to the agent in this scenario. We trained
the RL agent controlling the player for 3 million iterations, while monitoring its average
rewards over time. These iterations consisted of multiple games, called episodes. In each
episode, the agent and ﬂag are spawned at uniformly random locations in the ﬁeld. An
episode is completed when the player either successfully captures the ﬂag, or 1000 iterations
have passed. We selected the number of training iterations by examining when the average
reward per episode leveled near the maximum of 1000. The resulting trained agent was used
as a base model for proceeding teams that have a single player.

98

With a baseline agent trained, we added complexity to the task by introducing a second
team, as well as partitioning the ﬁeld into sides. In this new scenario, all of the rewards
in Table 7.1 were available. Each episode the side the agent defended switched to prevent
the agent from developing a preference to move the player in one direction. Training in
this scenario employed self-play, allowing the RL agent to learn against a previous version of
itself. That is, the new team was controlled by a static version of the baseline agent; however,
the RL agent continued to train. The RL agent was trained for an additional 7 millions steps
for a total of 10 million training steps. This process could be repeated, where the static team
is the one that was trained for 10 million steps, and the learning team continues to train
against its improved adversary. However, we only examined one iteration of self-play here.
The ﬁnal scenario we examined was the emergence of deterrence in the game. This sce-
nario is important for understanding how agents adapt to high-risk, high-reward situations,
aligning with our broader objective of examining strategic decision-making in complex en-
vironments. To this end, we improved the defending ability of the static agent by having
it following two rules: it would move directly towards the enemy at full speed if the enemy
was on its side of the ﬁeld, and otherwise, it would move randomly. To investigate the
threshold at which the RL agent would refrain from aggressive play, it was trained against
this improved defender, with a varying negative reward for being captured ranging from 0
to -10,000 points. The RL agent was trained for an additional 1 million iterations for each
value. As in previous scenarios, we ensured randomness in the initial placement of entities,
with each starting on their respective sides of the ﬁeld.

7.4 Results
Recall the simplest scenario we examined, where an RL agent was trained to control a single-
player team to capture a ﬂag anywhere in the ﬁeld. In Fig. 7.2, we show the progression of
average rewards per episode against training iterations in blue. Initially, the agent quickly
learns basic navigation, as evidenced by balancing negative and positive rewards within
approximately 500, 000 iterations. Within the subsequent 1.5 million iterations, the agent
reliably develops a strategy to move quickly towards the ﬂag. Given the random initial con-
ditions and step penalties, a perfect 1000 point reward is unattainable, leading to variation in
the maximum achievable reward. Nevertheless, the agent consistently attains rewards near
the theoretical maximum of 1000 points. Training was extended for an additional 1.5 million
iterations to reinforce the acquired behaviors and expose the agent to rarer edge cases. This
phase of training demonstrates the agent’s proﬁciency in mastering basic strategies, setting
a foundation for more complex scenarios.

Building on the understanding of the RL agent’s learning process, we next explore the

99

Figure 7.2: Average reward per episode versus training iteration of an RL agent learning to
capture a ﬂag with a single player. The RL agent learns basic navigation within 500,000
iterations and is able to balance negative and positive rewards. Within an additional 1
million iterations, the RL is reliably capturing the ﬂag.

speciﬁc actions undertaken by the trained agent in varying scenarios. Figure 7.3 illustrates
the agent’s action in four diﬀerent placements of the enemy ﬂag. Similarly, Fig. 7.4 demon-
strates the agent’s behavior in the absence of an enemy ﬂag. In these ﬁgures, the agent’s
actions are represented by arrows indicating the direction a player movement at each location,
with the enemy ﬂag marked as a red circle. In all cases, the agent has learned to prevent the
player from hitting the black boundary, as indicated by all actions by the boundary leading
the player away. In the scenario with no ﬂag, the agent’s general search strategy is evident,
characterized by moving the player in a clockwise pattern around the ﬁeld, spiralling around
central point. These behaviors highlight the agent’s adaptive strategies in diﬀerent game
environments.

In the next scenario, we shifted to training the agent against a static version of itself.
Figure 7.5 shows the actions the agent would supply a blue player at various ﬁeld positions,
represented by arrows. The red player, depicted as a smaller red circle aims to defend the
red ﬂag, while the RL agent is tasked with defending the blue ﬂag, both indicated by larger
colored circles. It is important to note that in this ﬁgure, we assume that both the red and
blue agents are stationary when determining these actions. However, in a live game scenario,
these agents likely to be in motion. Consequently, these predicted actions occur within a
parameter space the RL agent may not have previously encountered. This aspect should be
considered when interpreting the actions.

100

(a) Enemy ﬂag on the left half of ﬁeld.

(b) Enemy ﬂag on the top half of ﬁeld.

(c) Enemy ﬂag on the right half of ﬁeld.

(d) Enemy ﬂag on the bottom half of ﬁeld.

Figure 7.3: Actions of single player attacking a ﬂag. In each subﬁgure the arrows denote
the direction the attacking player would move. The black square is the boundary for the
ﬁeld, and the red circle is the enemy’s ﬂag. The agent has learned to avoid the boundary, as
denoted by all of the arrows pointing away from the boundary. Additionally, the agent has
generally learned how to move towards the ﬂag in various locations in the ﬁeld.

In Fig. 7.5a, the red player is positioned defensively. When the RL is on the oﬀensive, it
generally advances towards the enemy ﬂag. Its actions are unaﬀected by the defending red

101

Figure 7.4: Actions of a single player without a ﬂag. Similar to Fig. 7.3, the arrows denote
the direction the attacking player would move. The black square is the boundary for the
ﬁeld, and the red circle is the enemy’s ﬂag. When the agent is unable to detect a ﬂag, it
searches by moving clockwise around the ﬁeld around a central point.

player when distant from it. Interestingly, when the agent is near both the red player and
the center of the ﬁeld, it shifts towards its own ﬂag, indicating a defensive strategy. However,
in its own defensive position, the agent consistently moves towards its ﬂag. Notably, near
the bottom of the ﬁeld, the agent does not shift to an oﬀensive strategy, possibly indicating
a gap in its training.

Contrastingly, Fig. 7.5b places the red player in an oﬀensive position. In this scenario,
when the RL agent is on the oﬀensive it predominantly heads towards the enemy ﬂag, except
in areas adjacent to the red player near the center line. In a defensive position, the agent has
learned to target the attacking red player. This is evidence by arrows near the red player
directing towards it. These scenarios demonstrate the agent’s capacity for dynamic strategy
adjustments in response to diﬀerent threats.

In our ﬁnal analysis, we explored the emergence of deterrence by varying the negative
reward associated with the agent’s player being captured. Figure 7.6 presents the average
reward (depicted in blue) and average number of steps taken (in orange) across 30 simula-
tions, plotted against the varying penalties for being captured. As expected, with a minimal
penalty, the agent attempts to capture the ﬂag, indicated by the relatively low number of
steps coupled with a negative reward. However, as the penalty increases, a signiﬁcant shift
in the agent’s strategy is observed. Notably, at a penalty of approximately 8000 points, the

102

(a) If the red team is in a defensive position,
the blue player moves to capture the ﬂag when
they are in an oﬀensive position and defend their
ﬂag in a defensive one. There is a slight prefer-
ence for moving around the defensive red player.
Around the center line in front of the red player,
we see the agent’s actions leading to a defensive
position.

(b) If the red team is in an oﬀensive position,
the blue player moves to capture the agent when
they are in a defensive position, and move to
capture the ﬂag in an oﬀensive one. The blue
agent decides to move back to their side if they
are close to the center and near the red agent.

In these scenarios, the ﬁeld is
Figure 7.5: Actions of single player against single player.
partitioned into two sides by a grey dashed line, where the red team defends the left and
the blue team defends the right. The large circles denote the teams ﬂags, and the smaller
circle is a red player. The arrows denote the action a blue player would make, assuming all
entities are at rest.

agents average score raises to around -1000 points, with an average of 1000 steps taken. This
indicates a strategic change where the agent opts to remain on its side of the ﬁeld throughout
the game, eﬀectively exhibiting deterrence.

7.5 Conclusion
In this chapter, we successfully implemented a numerical version of capture the ﬂag and ap-
plied reinforcement learning to develop strategic gameplay. Our preliminary results demon-
strate the capability of training an RL agent to control a single player in progressively
complex scenarios.

Initially, the agent was tasked with capturing a randomly place ﬂag. Starting from this
basic task aided the RL agent in learning robust strategies.
In early attempts, an agent
always started on a speciﬁed side of the ﬁeld. Having never experienced the ﬂag behind it,

103

Figure 7.6: Average reward and number of steps versus penalty of being caught. The average
number of steps (orange) and reward (blue) for 30 simulations is shown versus the penalty
of being caught. When the penalty is low, the RL agent attempts to capture the enemy ﬂag.
However, once the penalty reaches approximately 8000, the agent is deterred from play and
spends the total 1000 timesteps on their side of the ﬁeld.

the agent learned to always move in a single direction regardless of its observations. This
insight inﬂuenced the initialization of subsequent scearios.

Once the baseline agent was trained to stay in bounds and capture a ﬂag, we moved to
training it against another player. When the new player followed a prescribed set of rules,
the RL excelled in adapting to oﬀensive rolls, however it struggled in subsequent defensive
roles. This challenge was partially mitigated through self-play that allowed the agent to
experience both oﬀensive and defensive roles during training and kept the opposing player
at a similar ability level. While this approach led to improved strategic behavior, such as
retreating to defend its ﬂag, it also highlighted areas for improvement. For example, as seen
in Fig. 7.5a, the agent occasionally fails to optimally avoid defenders.

The ﬁnal area we examined was the emergence of deterrence. By varying the negative
reward an RL agent was given for its player being captured, we were eventually able observe
it no longer play the game. This threshold was reached only at a relatively high penalty,
highlighting the agent’s initial propensity for risk-taking or goal-directed behavior despite
potential negative outcomes. This observation suggests that the agent’s decision making
process might require substantial negative reinforcement to override its objective-focused
strategies. Interestingly, it also implies that with further training iterations, a lower penalty
threshold might achieve similar deterrence eﬀects.

104

While the results we have been able to produce so far are informative, there are areas
that this research can be improved in the future. We sort these improvements into the
three broad categories of improved training, scenarios, and algorithms. Regarding train-
ing improvements, RL agents should be trained for many more iterations to provide ample
opportunities to develop and reﬁne sophisticated strategies, particularly through increased
self-play iterations. A signiﬁcant limitation of our work was the lack of hyperparameter
tuning, primarily driven by computational constraints. Nonetheless, future eﬀorts should
focus on exploring various neural network architectures, including adjustments in the num-
ber of hidden layers, their sizes, and activation functions. Additionally, ﬁne-tuning player
hyperparameters, such as the number and size of sensory rays, could lead to more nuanced
and eﬀective strategies.

In terms of scenario development, future research should explore crafting sculpted sce-
narios speciﬁcally designed to expedite the learning process of RL agents. In this work, the
starting locations of entities were randomized, which is eﬀective for broadly sample the state
space of the game. However, strategically selecting initial conditions to ensure diversity can
speed up the RL agent’s learning. This approach would expose the agent to a wide range
of situations in a more systematic manner, enhancing its ability to build strategies. Addi-
tionally, experimenting with varying team sizes and altering player attributes (e.g. player
hyperparmeters) presents an exciting opportunity. With these changes, we can dive deeper
into the dynamics of asymmetric conﬂict, examining how diﬀerences in team composition
and capabilities aﬀect strategic outcomes.

Finally, with respect to algorithms the more general multi-agent reinforcement learning
(MARL) should be explored. While we utilized single-agent reinforcement learning as a
starting point in this work, many complex real-world problems are not well represented in
the framework. Using MARL, both teams would be controlled by RL agents that learn simul-
taneously. Many RL algorithms can be generalized to multiple players. For example, [244]
examines extending trust region policy optimization, deep deterministic policy gradient, and
deep-q networks, while Kar et al. [245] discuss extending Q-learning to the multi-agent ver-
sion, QD-learning. Bu¸soniu et al. [246] provides a review many MARL algorithms. However,
transitioning from single-agent to multi-agent RL introduces unique challenges as outlined
by Nguyen el al. [247]. These areas include partial observability, non-stationarity, continu-
ous action spaces, multi-agent training schemes, and transfer learning. Our implementation
has already began to address partial partial observability and continuous spaces, providing
a solid foundation for tackling the other challenges in MARL.

In summary, this chapter lays the groundwork for addressing asymmetric conﬂict with
reinforcement learning through the use of a simulated capture the ﬂag game. Through many

105

scenarios, we have demonstrated the adaptability of RL agents in dynamic environments.
While are current results are promising, we have illuminated several areas for potential
improvement to further extend the capabilities of RL in simulations of conﬂict.

106

CHAPTER 8:

CONCLUSION
Throughout this dissertation, we explored areas of improvement for the ﬁeld of agent-based
modeling, namely data integration, policy evaluation, and the incorporation of machine
learning. These areas were examined through projects presented in pairs of chapters, where
each chapter consisted of a unique yet interconnected project.

Chapter 2 and 3 focused on using data to develop and reﬁne ABMs, speciﬁcally ap-
plied to the spread of chronic wasting disease in white tailed deer. These chapters serve as
comprehensive examples of data-driven model development and underscore the challenges in
modeling uncertainty. While our model produced simulations closely resembling real data,
the complexity of model creation and computational demands highlighted the need for more
eﬃcient methodologies.

The data integration projects underscored two major takeaways: the importance of ex-
ploratory data analysis (EDA) and sensitivity analysis in agent-based modeling. EDA was
instrumental in identifying subtle features present in our data that shaped the model’s de-
velopment. While integrating these features was challenging in our case, it may be more
straightforward in other applications. Future modeling eﬀorts should prioritize extensive
data examination to extract key features, thereby enhancing model accuracy across vari-
ous ﬁdelities. Similarly, sensitivity analysis (SA) plays an important role in model building
and assessment. SA aided in understanding the impact of parameters on model outcomes,
guiding decision on which features to include. For example, if a SA reveals that added
parameters have little to no eﬀect on important model outputs, you can safely omit those
parameters from the model. Moreover, SA can rank parameter importance, oﬀering insights
and understanding of the model that would otherwise be unattainable.

In Chapters 4 and 5, our research shifted towards employing ABMs for policy assessment.
Chapter 4 focused on evaluating strategies to mitigate the spread of disinformation across
social networks. We dedicated thousands of core hours to analyze and rank various strategies
across diﬀerent graph topologies and conducted a similar evaluation on a real social network.
However, the signiﬁcant computational demands posed challenges. These challenges led to
the focus of the subsequent chapter, where we explored ways to reduce the size of graphs while
preserving certain attributes. Such a process enhances the feasibility of applying models to
real-world social networks.

These projects highlighted the necessity of substantial computing resources for an in-
depth examination of policy implications using ABMs. Although our research demonstrated

107

techniques to lessen the demand for extensive resources, we were careful in our approach to
ensure these methods still accurately reﬂect real-world scenarios. This balance is crucial;
as we develop more eﬃcient computational strategies, it is imperative to maintain a high
degree of realism in our models. This ensures that the insights and recommendations derived
from ABMs remain relevant and applicable to actual policy-making. Future research should
continue to focus on optimizing computational eﬃciency while preserving the integrity and
realism of the models, thereby enabling a comprehensive and pragmatic evaluation of policies
in various domains.

Lastly, in Chapters 6 and 7, our focus shifted to the integration of machine learning within
ABMs, with a particular emphasis on models of conﬂict and deterrence. The ﬁrst chapter
in this segment laid the groundwork, enhancing existing mathematical models of conﬂict to
explicitly incorporate aspects of deterrence. Building on this, the following chapter presented
preliminary eﬀorts to apply reinforcement learning as the rules in agent-based models. This
approach was then used to examine the emergence of deterrence in a model of conﬂict, thus
illustrating how machine learning can contribute to the sophistication of ABMs in complex
scenarios.

While preliminary, these chapters illuminate the power of combining machine learning
with ABMs. This fusion not only enhances the analytical capabilities of ABMs but also
opens up new avenues for exploring intricate systems and behaviors. This advancement
points to a future where ABMs can be more eﬀectively used to simulate and understand
complex phenomena, oﬀering insights that might be diﬃcult to obtain through traditional
modeling methods. Future research can further reﬁne these methodologies, exploring their
applicability across a broader spectrum of disciplines and deepening our understanding of
the emergent properties within these systems.

The contributions of this dissertation to the ﬁeld of agent-based modeling highlight po-
tential avenues for future exploration. The dissertation demonstrates that integrating tra-
ditional modeling techniques with contemporary data and machine learning approaches can
lead to the development of more robust tools, better equipped to address the complexities
of today’s challenges.

108

BIBLIOGRAPHY

[1] Harry Jones. The recent large reduction in space launch cost. 48th International

Conference on Environmental Systems, 2018.

[2] Steven F. Railsback and Volker Grimm. Agent-Based and Individual-Based Modeling:

A Practical Introduction. Princeton University Press, 2011.

[3] C. M. Macal and M. J. North. Tutorial on agent-based modeling and simulation. In

Proceedings of the Winter Simulation Conference, 2005., pages 14 pp.–, 2005.

[4] C. M. Macal and M. J. North. Tutorial on agent-based modeling and simulation part 2:
How to model with agents. In Proceedings of the 2006 Winter Simulation Conference,
pages 73–83, 2006.

[5] John von Neumann and Arther Burks. Theory of Self-Replicating Automata. University

of Illinois Press, 1966.

[6] Martin Gardner. Mathematical Games The fantastic combinations of John Conway’s

new solitaire game “life”. Scientiﬁc American, pages 120–123, 1970.

[7] Stephen Wolfram. A New Kind of Science. Wolfram Media Inc., Champaign, Ilinois,

USA, 2002.

[8] Christopher Bone and Suzana Dragi´cevi´c. Simulation and validation of a reinforcement
learning agent-based model for multi-stakeholder forest management. Computers, En-
vironment and Urban Systems, 34(2):162–174, 2010.

[9] William Ogilvy Kermack, A. G. McKendrick, and Gilbert Thomas Walker. A con-
tribution to the mathematical theory of epidemics. Proceedings of the Royal Society
of London. Series A, Containing Papers of a Mathematical and Physical Character,
115(772):700–721, 1927.

[10] Thomas House, Geoﬀrey Davies, Leon Danon, and Matt J Keeling. A motif-based
approach to network epidemics. Bulletin of Mathematical Biology, 71:1693–1706, 2009.

[11] Jaewook Joo and Joel L Lebowitz. Pair approximation of the stochastic susceptible-
infected-recovered-susceptible epidemic model on the hypercubic lattice. Physical Re-
view E, 70(3):036114, 2004.

[12] Anirban Chakraborti, Ioane Muni Toke, Marco Patriarca, and Fr´ed´eric Abergel. Econo-
physics review: II. Agent-based models. Quantitative Finance, 11(7):1013–1041, 2011.

[13] M. Barbati, G. Bruno, and A. Genovese. Applications of agent-based models for opti-
mization problems: A literature review. Expert Systems with Applications, 39(5):6020
– 6028, 2012.

[14] Li An. Modeling human decisions in coupled human and natural systems: Review
of agent-based models. Ecological Modelling, 229:25 – 36, 2012. Modeling Human
Decisions.

109

[15] Federico Bianchi and Flaminio Squazzoni. Agent-based models in sociology. Wiley

Interdisciplinary Reviews: Computational Statistics, 7(4):284–306, 2015.

[16] Sean Barnes, Bruce Golden, and Stuart Price. Applications of agent-based model-
ing and simulation to healthcare operations management. In Handbook of healthcare
operations management: methods and applications, pages 45–74. Springer, 2013.

[17] Melissa Tracy, Magdalena Cerd´a, and Katherine M Keyes. Agent-based modeling in
public health: current applications and future directions. Annual review of public
health, 39:77–94, 2018.

[18] Charles M Macal. Everything you need to know about agent-based modelling and

simulation. Journal of Simulation, 10:144–156, 2016.

[19] Volker Grimm, Eloy Revilla, Uta Berger, Florian Jeltsch, Wolf M Mooij, Steven F
Railsback, Hans-Hermann Thulke, Jacob Weiner, Thorsten Wiegand, and Donald L
DeAngelis. Pattern-oriented modeling of agent-based complex systems: lessons from
ecology. science, 310(5750):987–991, 2005.

[20] Craig Loehle. A guide to increased creativity in research: inspiration or perspiration?

Bioscience, 40(2):123–129, 1990.

[21] Joshua M Epstein. Why model? Journal of artiﬁcial societies and social simulation,

11(4):12, 2008.

[22] Bruce Edmonds and Scott Moss. From kiss to kids–an ‘anti-simplistic’modelling ap-
proach. In International workshop on multi-agent systems and agent-based simulation,
pages 130–144. Springer, 2004.

[23] David O’Sullivan, Tom Evans, Steven Manson, Sara Metcalf, Arika Ligmann-Zielinska,
and Chris Bone. Strategic directions for agent-based modeling: avoiding the yaawn
syndrome. Journal of land use science, 11(2):177–187, 2016.

[24] William Rand and Roland T Rust. Agent-based modeling in marketing: Guidelines

for rigor. International Journal of research in Marketing, 28(3):181–193, 2011.

[25] William Rand. Theory-interpretable, data-driven agent-based modeling.

Social-

behavioral modeling for complex systems, pages 337–357, 2019.

[26] Srinivasan Venkatramanan, Bryan Lewis, Jiangzhuo Chen, Dave Higdon, Anil Vul-
likanti, and Madhav Marathe. Using data-driven agent-based models for forecasting
emerging infectious diseases. Epidemics, 22:43–49, 2018.

[27] Elizabeth Hunter, Brian Mac Namee, and John Kelleher. An open-data-driven agent-
based model to simulate infectious disease outbreaks. PloS one, 13(12):e0208775, 2018.

[28] Mazhar Sajjad, Karandeep Singh, Euihyun Paik, and Chang-Won Ahn. A data-driven
approach for agent-based modeling: Simulating the dynamics of family formation.
Journal of Artiﬁcial Societies and Social Simulation, 19(1):9, 2016.

110

[29] Haifeng Zhang, Yevgeniy Vorobeychik, Joshua Letchford, and Kiran Lakkaraju. Data-
driven agent-based modeling, with application to rooftop solar adoption. Autonomous
Agents and Multi-Agent Systems, 30:1023–1049, 2016.

[30] David J Butts, Noelle E Thompson, Sonja A Christensen, David M Williams, and
Michael S Murillo. Data-driven agent-based model building for animal movement
through exploratory data analysis. Ecological Modelling, 470:110001, 2022.

[31] David J Butts, Sam A Bollman, and Michael S Murillo. Mathematical modeling of
disinformation and eﬀectiveness of mitigation policies. Scientiﬁc Reports, 13(1):18735,
2023.

[32] Aldo Leopold. The Conservation Ethic. Journal of Forestry, 31(6):634–643, 10 1933.

[33] Hugh Dingle and V. Alistair Drake. What Is Migration? BioScience, 57(2):113–121,

02 2007.

[34] Navinder J. Singh, Andrew M. Allen, and G¨oran Ericsson. Quantifying migration
behaviour using net squared displacement approach: Clariﬁcations and caveats. PLOS
ONE, 11(3):1–20, 03 2016.

[35] Mevin B. Hooten, Henry R. Scharf, Trevor J. Heﬂey, Aaron T. Pearse, and Mitch D.
Weegman. Animal movement models for migratory individuals and groups. Methods
in Ecology and Evolution, 9(7):1692–1705, 2018.

[36] Sophie Bestley, Toby A. Patterson, Mark A. Hindell, and John S. Gunn. Feeding
ecology of wild migratory tunas revealed by archival tag records of visceral warming.
Journal of Animal Ecology, 77(6):1223–1233, 2008.

[37] B.F. Manly, L. McDonald, D.L. Thomas, T.L. McDonald, and W.P. Erickson. Re-
source Selection by Animals: Statistical Design and Analysis for Field Studies. Springer
Netherlands, 2007.

[38] Douglas H. Johnson. The comparison of usage and availability measurements for eval-

uating resource preference. Ecology, 61(1):65–71, 1980.

[39] Geert Aarts, Monique MacKenzie, Bernie McConnell, Mike Fedak, and Jason
Matthiopoulos. Estimating space-use and habitat preference from wildlife telemetry
data. Ecography, 31(1):140–160, 2008.

[40] D. John Anderson. The home range: A new nonparametric estimation technique.

Ecology, 63(1):103–112, 1982.

[41] D. B. Shepard, A. R. Kuhns, M. J. Dreslik, and C. A. Phillips. Roads as barriers
to animal movement in fragmented landscapes. Animal Conservation, 11(4):288–296,
2008.

[42] A. David M. Latham, M. Cecilia Latham, Mark S. Boyce, and Stan Boutin. Movement
responses by wolves to industrial linear features and their eﬀect on woodland caribou
in northeastern alberta. Ecological Applications, 21(8):2854–2865, 2011.

111

[43] Mevin B. Hooten, Devin S. Johnson, Brett T. McClintock, and Juan M. Morales.
Animal Movement Statistical Models For Telemetry Data. Taylor & Francis, 2017.

[44] Kim Whoriskey, Marie Auger-M`eth`e, Christoﬀer M. Albertsen, Frederick G. Who-
riskey, Thomas R. Binder, Charles C. Krueger, and Joanna Mills Flemming. A hidden
markov movement model for rapidly identifying behavioral states from animal tracks.
Ecology and Evolution, 7(7):2112–2121, 2017.

[45] Roland Langrock, Ruth King, Jason Matthiopoulos, Len Thomas, Daniel Fortin, and
Juan M. Morales. Flexible and practical modeling of animal telemetry data: hidden
markov models and extensions. Ecology, 93(11):2336–2342, 2012.

[46] Ann E. McKellar, Roland Langrock, Jeﬀrey R. Walters, and Dylan C. Kesler. Using
mixed hidden Markov models to examine behavioral states in a cooperatively breeding
bird. Behavioral Ecology, 26(1):148–157, 09 2014.

[47] Ian D. Jonsen, Ransom A. Myers, and Joanna Mills Flemming. META-ANALYSIS
OF ANIMAL MOVEMENT USING STATE-SPACE MODELS. Ecology, 84(11):3055–
3063, 2003.

[48] Ian D. Jonsen, Joanna Mills Flemming, and Ransom A. Myers. ROBUST STATE-
SPACE MODELING OF ANIMAL MOVEMENT DATA. Ecology, 86(11):2874–2880,
2005.

[49] Toby A. Patterson, Len Thomas, Chris Wilcox, Otso Ovaskainen, and Jason
Matthiopoulos. State-space models of individual animal movement. Trends in Ecology
& Evolution, 23(2):87 – 94, 2008.

[50] Haiganoush K. Preisler, Alan A. Ager, and Michael J. Wisdom. Analyzing animal

movement patterns using potential functions. Ecosphere, 4(3):art32, 2013.

[51] Th`eo Michelot, Pierre Gloaguen, Paul G. Blackwell, and Marie-Pierre ˜Atienne. The
langevin diﬀusion as a continuous-time model of animal movement and habitat selec-
tion. Methods in Ecology and Evolution, 10(11):1894–1907, 2019.

[52] John W Tukey et al. Exploratory data analysis, volume 2. Reading, Mass., 1977.

[53] Frederick Hartwig and Brian E Dearing. Exploratory data analysis. Sage, 1979.

[54] Matthew B Miles and A Michael Huberman. Qualitative data analysis: An expanded

sourcebook. sage, 1994.

[55] Andrew M. Edwards. Overturning conclusions of l`evy ﬂight movement patterns by

ﬁshing boats and foraging animals. Ecology, 92(6):1247–1257, 2011.

[56] Andrew M. Edwards, Richard A. Phillips, Nicholas W. Watkins, Mervyn P. Free-
man, Eugene J. Murphy, Vsevolod Afanasyev, Sergey V. Buldyrev, M. G. E. da Luz,
E. P. Raposo, H. Eugene Stanley, and Gandhimohan M. Viswanathan. Revisiting
l´evy ﬂight search patterns of wandering albatrosses, bumblebees and deer. Nature,
449(7165):1044–1048, 2007.

112

[57] Andrew M. Edwards. Using likelihood to test for l`evy ﬂight search patterns and for
general power-law distributions in nature. Journal of Animal Ecology, 77(6):1212–1222,
2008.

[58] Erhard Rahm and Hong Hai Do. Data cleaning: Problems and current approaches.

IEEE Data Eng. Bull., 23(4):3–13, 2000.

[59] Tamraparni Dasu and Theodore Johnson. Exploratory data mining and data cleaning,

volume 479. John Wiley & Sons, 2003.

[60] Jonathan Schwabish. Better Data Visualizations: A Guide for Scholars, Researchers,

and Wonks. Columbia University Press, 2021.

[61] Robert M Edsall. The parallel coordinate plot in action: design and use for geographic

visualization. Computational Statistics & Data Analysis, 43(4):605–619, 2003.

[62] Takayuki Itoh, Ashnil Kumar, Karsten Klein, and Jinman Kim. High-dimensional data
visualization by interactive construction of low-dimensional parallel coordinate plots.
Journal of Visual Languages & Computing, 43:1–13, 2017.

[63] Andrej Gisbrecht and Barbara Hammer. Data visualization by nonlinear dimensional-
ity reduction. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery,
5(2):51–73, 2015.

[64] Samuel Kaski and Jaakko Peltonen. Dimensionality reduction for data visualization

[applications corner]. IEEE signal processing magazine, 28(2):100–104, 2011.

[65] John D. Hunter. Matplotlib: A 2d graphics environment. Computing in Science

Engineering, 9(3):90–95, 2007.

[66] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blon-
del, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau,
M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python.
Journal of Machine Learning Research, 12:2825–2830, 2011.

[67] Michael L. Waskom. seaborn: statistical data visualization. Journal of Open Source

Software, 6(60):3021, 2021.

[68] Benjamin Bengfort and Rebecca Bilbro. Yellowbrick: Visualizing the scikit-learn model

selection process. Journal of Open Source Software, 4(35):1075, 2019.

[69] Amy Quinn. Inﬂuences of movement behavior and space use in evaluating disease risk
amoung white-tailed deer in central new york. Syracuse, New York: State University
of New York College of Environmental Science and Forestry, 2010.

[70] David M. Williams, Amy C. Dechen Quinn, and William F. Porter. Informing disease
models with temporal and spatial contact structure among gps-collared individuals in
wild populations. PLOS ONE, 9, 01 2014.

113

[71] Mark S Boyce, Pierre R Vernier, Scott E Nielsen, and Fiona K.A Schmiegelow. Eval-
uating resource selection functions. Ecological Modelling, 157(2):281–300, 2002.

[72] Toby A. Patterson, Marinelle Basson, Mark V. Bravington, and John S. Gunn. Classi-
fying movement behaviour in relation to environmental conditions using hidden markov
models. Journal of Animal Ecology, 78(6):1113–1123, 2009.

[73] D. R. Brillinger, H. K. Preisler, A. A. Ager, and J. G. Kie. The Use Of Potential
Functions In Modelling Animal Movement, pages 385–409. Springer New York, New
York, NY, 2012.

[74] David R. Brillinger, Haiganoush K. Preisler, Alan A. Ager, John G. Kie, and Brent S.
Stewart. Employing stochastic diﬀerential equations to model wildlife motion. Bulletin
of the Brazilian Mathematical Society, 33(3):385–408, 2002.

[75] A. Quinn Dechen, David Williams, W. Porter, M. Smith, F. DeSantis, and McNulty
F. Risk assessment of a chronic wasting disease outbreak in new york: Final report.
New York State Department of Environmental Conservation, 2009.

[76] Daniel Fortin, Hawthorne L. Beyer, Mark S. Boyce, Douglas W. Smith, Thierry Duch-
esne, and Julie S. Mao. Wolves inﬂuence elk movements: Behavior shapes a trophic
cascade in yellowstone national park. Ecology, 86(5):1320–1330, 2005.

[77] Robert Tibshirani, Michael Saunders, Saharon Rosset, Ji Zhu, and Keith Knight.
Sparsity and smoothness via the fused lasso. Journal of the Royal Statistical Society:
Series B (Statistical Methodology), 67(1):91–108, 2005.

[78] Steven Diamond and Stephen Boyd. CVXPY: A Python-embedded modeling language
for convex optimization. Journal of Machine Learning Research, 17(83):1–5, 2016.

[79] Akshay Agrawal, Robin Verschueren, Steven Diamond, and Stephen Boyd. A rewriting
system for convex optimization problems. Journal of Control and Decision, 5(1):42–60,
2018.

[80] Tak chung Fu. A review on time series data mining. Engineering Applications of

Artiﬁcial Intelligence, 24(1):164 – 181, 2011.

[81] Christos Berberidis, Walid G. Aref, Mikhail Atallah, Ioannis Vlahavas, and Ahmed K.
Elmagarmid. Multiple and parital periodicity in mining in time series databases. In
Proceedings of the 15th European Conerene on Artiﬁcial Intelligence, 2002.

[82] St´ephane Dray, Manuela Royer-Carenzi, and Cl´ement Calenge. The exploratory analy-
sis of autocorrelation in animal-movement studies. Ecological Research, 25(3):673–681,
2010.

[83] Guillaume P´eron, Chris H. Fleming, Rogerio C. de Paula, and Justin M. Calabrese.
Uncovering periodic patterns of space use in animal tracking data with periodograms,
including a new algorithm for the lomb-scargle periodogram and improved randomiza-
tion tests. Movement Ecology, 4(1):19, 2016.

114

[84] Jacob T. VanderPlas. Understanding the lomb–scargle periodogram. The Astrophysical

Journal Supplement Series, 236(1):16, may 2018.

[85] Devin S. Johnson, Joshua M. London, Mary-Anne Lea, and John W. Durban.
Continuous-time correlated random walk model for animal telemetry data. Ecology,
89(5):1208–1215, 2008.

[86] Tsutomu T. Takeuchi. Constructing a bivariate distribution function with given
marginals and correlation: application to the galaxy luminosity function. Monthly
Notices of the Royal Astronomical Society, 406(3):1830–1840, 08 2010.

[87] Pravin K. Trivedi and David M. Zimmer. Copula Modeling: An Introduction for

Practitioners. Foundations and Trends(R) in Econometrics, 1(1):1–111, April 2007.

[88] Yajie Zou, Xinzhi Zhong, Jinjun Tang, Xin Ye, Lingtao Wu, Muhammad Ijaz, and
Yinhai Wang. A copula-based approach for accommodating the underreporting eﬀect
in wildlife-vehicle crash analysis. Sustainability, 11(2), 2019.

[89] Justin T. French, Hsiao-Hsuan Wang, William E. Grant, and John M. Tomeˇcek. Dy-
namics of animal joint space use: a novel application of a time series approach. Move-
ment Ecology, 7(1):38, 2019.

[90] Marti J. Anderson, Perry de Valpine, Andrew Punnett, and Arden E. Miller. A path-
way for multivariate analysis of ecological communities using copulas. Ecology and
Evolution, 9(6):3276–3294, 2019.

[91] B. J. Worton. Kernel methods for estimating the utilization distribution in home-range

studies. Ecology, 70(1):164–168, 1989.

[92] D. Erran Seaman and Roger A. Powell. An evaluation of the accuracy of kernel density

estimators for home range analysis. Ecology, 77(7):2075–2085, 1996.

[93] G. E. Uhlenbeck and L. S. Ornstein. On the theory of the brownian motion. Phys.

Rev., 36:823–841, Sep 1930.

[94] Volker Grimm, Steven F Railsback, Christian E Vincenot, Uta Berger, Cara Gallagher,
Donald L DeAngelis, Bruce Edmonds, Jiaqi Ge, Jarl Giske, Juergen Groeneveld, et al.
The odd protocol for describing agent-based and other simulation models: A second
update to improve clarity, replication, and structural realism. Journal of Artiﬁcial
Societies and Social Simulation, 23(2), 2020.

[95] Noelle E Thompson, David J Butts, Michael S Murillo, David M Williams, Sonja A
Christensen, William F Porter, and Gary J Roloﬀ. An individual-based model for
direct and indirect transmission of chronic wasting disease in free-rangin white-tailed
deer. Ecological Modelling, (under review).

[96] Ryan G McClarren, Penrose McClarren, and RL Penrose. Uncertainty quantiﬁcation

and predictive computational science. Springer, 2018.

115

[97] Arika Ligmann-Zielinska, Peer-Olaf Siebers, Nicholas Magliocca, Dawn C Parker,
Volker Grimm, Jing Du, Martin Cenek, Viktoriia Radchuk, Nazia N Arbab, Sheng
Li, et al. ‘one size does not ﬁt all’: a roadmap of purpose-driven mixed-method path-
ways for sensitivity analysis of agent-based models. Journal of Artiﬁcial Societies and
Social Simulation, 23(1), 2020.

[98] Andrea Saltelli, Stefano Tarantola, and KP-S Chan. A quantitative model-independent
method for global sensitivity analysis of model output. Technometrics, 41(1):39–56,
1999.

[99] J´erˆome Morio. Global and local sensitivity analysis methods for a physical system.

European journal of physics, 32(6):1577, 2011.

[100] Arika Ligmann-Zielinska.

Spatially-explicit sensitivity analysis of an agent-based
model of land use change. International Journal of Geographical Information Science,
27(9):1764–1781, 2013.

[101] Bertrand Iooss and Paul Lemaˆıtre. A review on global sensitivity analysis methods.
Uncertainty management in simulation-optimization of complex systems: algorithms
and applications, pages 101–122, 2015.

[102] Andrea Saltelli, Stefano Tarantola, Francesca Campolongo, Marco Ratto, et al. Sen-
sitivity analysis in practice: a guide to assessing scientiﬁc models, volume 1. Wiley
Online Library, 2004.

[103] Emanuele Borgonovo and Elmar Plischke. Sensitivity analysis: A review of recent

advances. European Journal of Operational Research, 248(3):869–887, 2016.

[104] Andrea Saltelli, Paola Annoni, Ivano Azzini, Francesca Campolongo, Marco Ratto, and
Stefano Tarantola. Variance based sensitivity analysis of model output. design and es-
timator for the total sensitivity index. Computer physics communications, 181(2):259–
270, 2010.

[105] IM Sobo´l. Sensitivity estimates for nonlinear mathematical models. Math. Model.

Comput. Exp., 1:407, 1993.

[106] Toshimitsu Homma and Andrea Saltelli.

Importance measures in global sensitivity
analysis of nonlinear models. Reliability Engineering & System Safety, 52(1):1–17,
1996.

[107] Max D Morris. Factorial sampling plans for preliminary computational experiments.

Technometrics, 33(2):161–174, 1991.

[108] Adam Badawy, Emilio Ferrara, and Kristina Lerman. Analyzing the digital traces
of political manipulation: The 2016 russian interference twitter campaign.
In 2018
IEEE/ACM International Conference on Advances in Social Networks Analysis and
Mining (ASONAM), pages 258–265, 2018.

116

[109] Adam Fourney, Miklos Z Racz, Gireeja Ranade, Markus Mobius, and Eric Horvitz. Ge-
ographic and temporal trends in fake news consumption during the 2016 us presidential
election. In CIKM, volume 17, pages 6–10, 2017.

[110] Sahil Loomba, Alexandre de Figueiredo, Simon J Piatek, Kristen de Graaf, and Heidi J
Larson. Measuring the impact of covid-19 vaccine misinformation on vaccination intent
in the uk and usa. Nature human behaviour, 5(3):337–348, 2021.

[111] Talha Burki. Vaccine misinformation and social media. The Lancet Digital Health,

1(6):e258–e259, 2019.

[112] Warren Cornwall. Oﬃcials gird for a war on vaccine misinformation, 2020.

[113] Matthew W. Hughey. The who and why of qanon’s rapid rise. New Labor Forum,

30(3):76–87, 2021.

[114] Derek du Preez. Chatgpt has the potential to spread misinformation ‘at unprece-
dented scale’. https://diginomica.com/chatgpt-has-potential-spread-misinformation-
unprecedented-scale, 2023.

[115] Jack Brewster, Lorenzo Arvanitis, and McKenzie Sadeghi. Misinformation monitor:
January 2023. https://www.newsguardtech.com/misinformation-monitor/jan-2023/,
2023.

[116] Chatbots trigger next misinformation nightmare. https://www.axios.com/2023/02/21

/chatbots-misinformation-nightmare-chatgpt-ai, 2023.

[117] Josh A. Goldstein, Girish Sastry, Micah Musser, Renee DiResta, Matthew Gentzel, and
Katerina Sedova. Generative language models and automated inﬂuence operations:
Emerging threats and potential mitigations, 2023.

[118] Monther Aldwairi and Ali Alwahedi. Detecting fake news in social media networks.
Procedia Computer Science, 141:215–222, 2018. The 9th International Conference on
Emerging Ubiquitous Systems and Pervasive Networks (EUSPN-2018) / The 8th Inter-
national Conference on Current and Future Trends of Information and Communication
Technologies in Healthcare (ICTH-2018) / Aﬃliated Workshops.

[119] Supanya Aphiwongsophon and Prabhas Chongstitvatana. Detecting fake news with
In 2018 15th International Conference on Electrical En-
machine learning method.
gineering/Electronics, Computer, Telecommunications and Information Technology
(ECTI-CON), pages 528–531, 2018.

[120] Jun Lin, Glenna Tremblay-Taylor, Guanyi Mou, Di You, and Kyumin Lee. Detecting
fake news articles. In 2019 IEEE International Conference on Big Data (Big Data),
pages 3021–3025, 2019.

[121] Stephanie Preston, Anthony Anderson, David J. Robertson, Mark P. Shephard, and
Narisong Huhe. Detecting fake news on facebook: The role of emotional intelligence.
PLOS ONE, 16(3):1–13, 03 2021.

117

[122] Kai Shu, Amrita Bhattacharjee, Faisal Alatawi, Tahora H Nazer, Kaize Ding, Man-
sooreh Karami, and Huan Liu. Combating disinformation in a social media age. Wiley
Interdisciplinary Reviews: Data Mining and Knowledge Discovery, 10(6):e1385, 2020.

[123] Florian Saurwein and Charlotte Spencer-Smith. Combating disinformation on social
media: Multilevel governance and distributed accountability in europe. Digital Jour-
nalism, 8(6):820–841, 2020.

[124] Joanna M Burkhardt. Combating fake news in the digital age, volume 53. American

Library Association Chicago, IL, 2017.

[125] Karishma Sharma, Feng Qian, He Jiang, Natali Ruchansky, Ming Zhang, and Yan
Liu. Combating fake news: A survey on identiﬁcation and mitigation techniques.
ACM Transactions on Intelligent Systems and Technology (TIST), 10(3):1–42, 2019.

[126] Tools that ﬁght disinformation online. https://www.rand.org/research/projects/truth-

decay/ﬁghting-disinformation/search.html.

[127] JA Gallo and CY Cho. Social media: Misinformation and content moderation issues

for congress. Congressional Research Service Report, 46662, 2021.

[128] David MJ Lazer, Matthew A Baum, Yochai Benkler, Adam J Berinsky, Kelly M Green-
hill, Filippo Menczer, Miriam J Metzger, Brendan Nyhan, Gordon Pennycook, David
Rothschild, et al. The science of fake news. Science, 359(6380):1094–1096, 2018.

[129] Antino Kim, Patricia L Moravec, and Alan R Dennis. Combating fake news on social
media with source ratings: The eﬀects of user and expert reputation ratings. Journal
of Management Information Systems, 36(3):931–968, 2019.

[130] Ceren Budak, Divyakant Agrawal, and Amr El Abbadi. Limiting the spread of mis-
information in social networks. In Proceedings of the 20th international conference on
World wide web, pages 665–674, 2011.

[131] Xinran He, Guojie Song, Wei Chen, and Qingye Jiang. Inﬂuence blocking maximization
in social networks under the competitive linear threshold model. In Proceedings of the
2012 siam international conference on data mining, pages 463–474. SIAM, 2012.

[132] Mei Li, Xiang Wang, Kai Gao, and Shanshan Zhang. A survey on information diﬀusion

in online social networks: Models and methods. Information, 8(4), 2017.

[133] David Kempe, Jon Kleinberg, and ´Eva Tardos. Maximizing the spread of inﬂuence
through a social network. In Proceedings of the ninth ACM SIGKDD international
conference on Knowledge discovery and data mining, pages 137–146, 2003.

[134] Thomas M Liggett et al. Stochastic interacting systems: contact, voter and exclusion

processes, volume 324. springer science & Business Media, 1999.

[135] Giordano De Marzo, Andrea Zaccaria, and Claudio Castellano. Emergence of polar-
ization in a voter model with personalized information. Physical Review Research,
2(4):043117, 2020.

118

[136] Lu´ıs Carlos F Latoski, WG Dantas, and Jeferson J Arenzon. Curvature-driven growth
and interfacial noise in the voter model with self-induced zealots. Physical Review E,
106(1):014121, 2022.

[137] Robert Axelrod. The dissemination of culture: A model with local convergence and

global polarization. Journal of conﬂict resolution, 41(2):203–226, 1997.

[138] William Ogilvy Kermack and Anderson G McKendrick. A contribution to the math-
ematical theory of epidemics. Proceedings of the royal society of london. Series A,
Containing papers of a mathematical and physical character, 115(772):700–721, 1927.

[139] Daryl J Daley and David G Kendall. Epidemics and rumours. Nature, 204(4963):1118–

1118, 1964.

[140] Daniel P Maki and Maynard Thompson. Mathematical models and applications: with

emphasis on the social, life, and management sciences. Prentice Hall, 1973.

[141] Jos´e RC Piqueira, Mauro Zilbovicius, and Cristiane M Batistela. Daley–kendal mod-
els in fake-news scenario. Physica A: Statistical Mechanics and its Applications,
548:123406, 2020.

[142] Saeed Jamalzadeh, Kash Barker, Andr´es D Gonz´alez, and Sridhar Radhakrishnan.
Protecting infrastructure performance from disinformation attacks. Scientiﬁc Reports,
12(1):12707, 2022.

[143] Robert Axelrod, Joshua J Daymude, and Stephanie Forrest. Preventing extreme polar-
ization of political attitudes. Proceedings of the National Academy of Sciences, 118(50),
2021.

[144] Luc Steels. A self-organizing spatial vocabulary. Artiﬁcial life, 2(3):319–332, 1995.

[145] Xiang Niu, Casey Doyle, Gyorgy Korniss, and Boleslaw K Szymanski. The impact of
variable commitment in the naming game on consensus formation. Scientiﬁc reports,
7(1):1–11, 2017.

[146] Jierui Xie, Sameet Sreenivasan, Gyorgy Korniss, Weituo Zhang, Chjan Lim, and
Boleslaw K Szymanski. Social consensus through the inﬂuence of committed minorities.
Physical Review E, 84(1):011130, 2011.

[147] Jierui Xie, Jeﬀrey Emenheiser, Matthew Kirby, Sameet Sreenivasan, Boleslaw K Szy-
manski, and Gyorgy Korniss. Evolution of opinions on social networks in the presence
of competing committed groups. PLoS One, 7(3):e33215, 2012.

[148] Damon Centola, Joshua Becker, Devon Brackbill, and Andrea Baronchelli. Experi-
mental evidence for tipping points in social convention. Science, 360(6393):1116–1119,
jun 2018.

[149] Sreeja Nair, Kin Wai Ng, Adriana Iamnitchi, and John Skvoretz. Diﬀusion of social
conventions across polarized communities: an empirical study. Social Network Analysis
and Mining, 11:1–17, 2021.

119

[150] William W Hackborn, Tetiana Reznychenko, and Yihang Zhang. Consensus building

by committed agents. CODEE Journal, 12(1):2, 2019.

[151] David Galehouse, Tommy Nguyen, Sameet Sreenivasan, Omar Lizardo, G Korniss,
and B Szymanski. Impact of network connectivity and agent commitment on spread
of opinions in social networks. In Proceedings of the 5th International Conference on
Applied Human Factors and Ergonomics, pages 2318–2329, 2014.

[152] Andrew M Thompson, Boleslaw K Szymanski, and Chjan C Lim. Propensity and
stickiness in the naming game: Tipping fractions of minorities. Physical review E,
90(4):042809, 2014.

[153] Casey Doyle, Sameet Sreenivasan, Boleslaw K Szymanski, and Gyorgy Korniss. Social
consensus and tipping points with opinion inertia. Physica A: Statistical Mechanics
and its Applications, 443:316–323, 2016.

[154] Mauro Mobilia. Commitment versus persuasion in the three-party constrained voter

model. Journal of Statistical Physics, 151(1):69–91, 2013.

[155] Weituo Zhang, Chjan Lim, and Boleslaw K Szymanski. Analytic treatment of tipping
points for social consensus in large random networks. Physical Review E, 86(6):061134,
2012.

[156] Dina Mistry, Qian Zhang, Nicola Perra, and Andrea Baronchelli. Committed activists
and the reshaping of status-quo social consensus. Physical Review E, 92(4):042805,
2015.

[157] Benedek Rozemberczki and Rik Sarkar. Characteristic Functions on Graphs: Birds
of a Feather, from Statistical Descriptors to Parametric Models.
In Proceedings of
the 29th ACM International Conference on Information and Knowledge Management
(CIKM ’20), page 1325–1334. ACM, 2020.

[158] Aric Hagberg, Pieter Swart, and Daniel S Chult. Exploring network structure, dynam-
ics, and function using networkx. Technical report, Los Alamos National Lab.(LANL),
Los Alamos, NM (United States), 2008.

[159] Yifan Hu. Eﬃcient, high-quality force-directed graph drawing. Mathematica journal,

10(1):37–71, 2005.

[160] How to

tailor

covid-19

vaccine

information

to

your

speciﬁc

audience.

https://www.cdc.gov/vaccines/covid-19/hcp/tailoring-information.html.

[161] Mark Newman. Networks. Oxford university press, 2018.

[162] Bert Vogelstein, David Lane, and Arnold J. Levine. Surﬁng the p53 network. Na-
ture, 408(6810):307–310, November 2000. Number: 6810 Publisher: Nature Publishing
Group.

120

[163] Jingchun Chen and Bo Yuan. Detecting functional modules in the yeast protein–protein

interaction network. Bioinformatics, 22(18):2283–2290, September 2006.

[164] Albert-Laszlo Barabasi and Zoltan N Oltvai. Network biology: understanding the cell’s

functional organization. Nature reviews genetics, 5(2):101–113, 2004.

[165] Matthew O. Jackson. Social and economic networks. Princeton University Press,

Princeton, N.J. Woodstock, 2011.

[166] Morris H DeGroot. Reaching a consensus. Journal of the American Statistical associ-

ation, 69(345):118–121, 1974.

[167] Nicholas Economides. The economics of networks. International Journal of Industrial

Organization, 14(6):673–699, October 1996.

[168] W Brian Arthur. Increasing returns and path dependence in the economy. University

of michigan Press, 1994.

[169] Romualdo Pastor-Satorras and Alessandro Vespignani. Epidemic spreading in scale-

free networks. Physical review letters, 86(14):3200, 2001.

[170] Istv´an Z Kiss, Joel C Miller, P´eter L Simon, et al. Mathematics of epidemics on

networks. Cham: Springer, 598:31, 2017.

[171] Thomas House. Modelling epidemics on networks. Contemporary Physics, 53(3):213–

225, 2012.

[172] Sebastiano A Delre, Wander Jager, Tammo HA Bijmolt, and Marco A Janssen. Will
it spread or not? the eﬀects of social inﬂuences and network topology on innovation
diﬀusion. Journal of Product Innovation Management, 27(2):267–282, 2010.

[173] Daron Acemoglu, Asuman Ozdaglar, and Ali ParandehGheibi. Spread of (mis) infor-
mation in social networks. Games and Economic Behavior, 70(2):194–227, 2010.

[174] Ryan A. Rossi and Nesreen K. Ahmed. The network data repository with interactive

graph analytics and visualization. In AAAI, 2015.

[175] Jure Leskovec and Andrej Krevl. SNAP Datasets: Stanford large network dataset

collection. http://snap.stanford.edu/data, June 2014.

[176] Marinka Zitnik, Rok Sosiˇc, Sagar Maheshwari, and Jure Leskovec. BioSNAP Datasets:
Stanford biomedical network dataset collection. http://snap.stanford.edu/biodata, Au-
gust 2018.

[177] Srijan Sengupta. Statistical Network Analysis: Past, Present, and Future, October

2023. arXiv:2311.00122 [stat].

121

[178] Grzegorz Malewicz, Matthew H. Austern, Aart J.C Bik, James C. Dehnert, Ilan Horn,
Naty Leiser, and Grzegorz Czajkowski. Pregel: a system for large-scale graph process-
ing. In Proceedings of the 2010 ACM SIGMOD International Conference on Manage-
ment of data, SIGMOD ’10, pages 135–146, New York, NY, USA, June 2010. Associ-
ation for Computing Machinery.

[179] Robert Ryan McCune, Tim Weninger, and Greg Madey. Thinking Like a Vertex: A
Survey of Vertex-Centric Frameworks for Large-Scale Distributed Graph Processing.
ACM Computing Surveys, 48(2):25:1–25:39, October 2015.

[180] Nadathur Satish, Changkyu Kim, Jatin Chhugani, Hideki Saito, Rakesh Krishnaiyer,
Mikhail Smelyanskiy, Milind Girkar, and Pradeep Dubey. Can traditional program-
ming bridge the Ninja performance gap for parallel computing applications? ACM
SIGARCH Computer Architecture News, 40(3):440–451, June 2012.

[181] Nadathur Satish, Narayanan Sundaram, Md. Mostofa Ali Patwary, Jiwon Seo, Jongsoo
Park, M. Amber Hassaan, Shubho Sengupta, Zhaoming Yin, and Pradeep Dubey.
Navigating the maze of graph analytics frameworks using massive graph datasets. In
Proceedings of the 2014 ACM SIGMOD International Conference on Management of
Data, SIGMOD ’14, pages 979–990, New York, NY, USA, June 2014. Association for
Computing Machinery.

[182] Paul Erd˝os, Alfr´ed R´enyi, et al. On the evolution of random graphs. Publ. math. inst.

hung. acad. sci, 5(1):17–60, 1960.

[183] Albert-L´aszl´o Barab´asi and R´eka Albert. Emergence of scaling in random networks.

science, 286(5439):509–512, 1999.

[184] David J Earl and Michael W Deem. Parallel tempering: Theory, applications, and
new perspectives. Physical Chemistry Chemical Physics, 7(23):3910–3916, 2005.

[185] Robert H Swendsen and Jian-Sheng Wang. Replica monte carlo simulation of spin-

glasses. Physical review letters, 57(21):2607, 1986.

[186] Gareth O Roberts and Jeﬀrey S Rosenthal. Optimal scaling for various metropolis-

hastings algorithms. Statistical science, 16(4):351–367, 2001.

[187] Andrew Gelman, Walter R Gilks, and Gareth O Roberts. Weak convergence and op-
timal scaling of random walk metropolis algorithms. The annals of applied probability,
7(1):110–120, 1997.

[188] Hamsterster. Hamsterster social network. http://www.hamsterster.com.

[189] Mathieu Bastian, Sebastien Heymann, and Mathieu Jacomy. Gephi: an open source
software for exploring and manipulating networks. In Third international AAAI con-
ference on weblogs and social media, 2009.

122

[190] Yves F Atchad´e, Gareth O Roberts, and Jeﬀrey S Rosenthal. Towards optimal scaling
of metropolis-coupled markov chain monte carlo. Statistics and Computing, 21:555–
568, 2011.

[191] Paul K Huth. Deterrence and international conﬂict: Empirical ﬁndings and theoretical

debates. Annual Review of Political Science, 2(1):25–48, 1999.

[192] Moshe Kress. Modeling armed conﬂicts. Science, 336(6083):865–869, 2012.

[193] Gerardo Minguela-Castro, Ruben Heradio, and Carlos Cerrada. Automated support

for battle decision making. Military Operations Research, 27(4):5–24, 2022.

[194] Thomas C Schelling. Arms and inﬂuence. In Strategic Studies, pages 96–114. Rout-

ledge, 2008.

[195] Frederick William Lanchester. Aircraft in warfare: The dawn of the fourth arm. Con-

stable limited, 1916.

[196] Manvi Sahni and Sumanta Kumar Das. Performance of maximum likelihood estimator
for ﬁtting lanchester equations on kursk battle data. Journal of Battleﬁeld Technology,
18(2):23–30, 2015.

[197] Thomas W Lucas and Turker Turkes. Fitting lanchester equations to the battles of

kursk and ardennes. Naval Research Logistics (NRL), 51(1):95–116, 2004.

[198] Ian R Johnson and Niall J MacKay. Lanchester models and the battle of britain. Naval

Research Logistics (NRL), 58(3):210–222, 2011.

[199] Jerome Bracken. Lanchester models of the ardennes campaign. Naval Research Logis-

tics (NRL), 42(4):559–577, 1995.

[200] Moshe Kress, Jonathan P. Caulkins, Gustav Feichtinger, Dieter Grass, and Andrea
Seidl. Lanchester model for three-way combat. European Journal of Operational Re-
search, 264(1):46–54, 2018.

[201] Moshe Kress. Lanchester models for irregular warfare. Mathematics, 8(5):737, 2020.

[202] Eduardo Gonz´alez and Marcelo Villena. Spatial lanchester models. European Journal

of Operational Research, 210(3):706–715, 2011.

[203] William R Caspary. Richardson’s model of arms races: description, critique, and an

alternative model. International Studies Quarterly, 11(1):63–88, 1967.

[204] Norman Z Alcock and Keith Lowe. The vietnam war as a richardson process. Journal

of Peace Research, 6(2):105–111, 1969.

[205] Jean-Christian Lambelet. A dynamic model of the arms race in the Middle East,

1953-1965. Society for General Systems Research, 1971.

123

[206] Seymour J Deitchman. A lanchester model of guerrilla warfare. Operations Research,

10(6):818–827, 1962.

[207] Hannes Risken and Hannes Risken. Fokker-planck equation. Springer, 1996.

[208] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy,
David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan
Bright, St´efan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman,
Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J
Carey, ˙Ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef
Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M.
Archibald, Antˆonio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0
Contributors. SciPy 1.0: Fundamental Algorithms for Scientiﬁc Computing in Python.
Nature Methods, 17:261–272, 2020.

[209] Takuya Iwanaga, William Usher, and Jonathan Herman. Toward SALib 2.0: Ad-
vancing the accessibility and interpretability of global sensitivity analyses. Socio-
Environmental Systems Modelling, 4:18155, May 2022.

[210] Jon Herman and Will Usher. SALib: An open-source python library for sensitivity

analysis. The Journal of Open Source Software, 2(9), jan 2017.

[211] Francesca Pianosi and Thorsten Wagener. A simple and eﬃcient method for global sen-
sitivity analysis based on cumulative distribution functions. Environmental Modelling
& Software, 67:1–11, 2015.

[212] Mark Towers, Jordan K. Terry, Ariel Kwiatkowski, John U. Balis, Gianluca de Cola,
Tristan Deleu, Manuel Goul˜ao, Andreas Kallinteris, Arjun KG, Markus Krimmel, Ro-
drigo Perez-Vicente, Andrea Pierr´e, Sander Schulhoﬀ, Jun Jet Tai, Andrew Tan Jin
Shen, and Omar G. Younis. Gymnasium, March 2023.

[213] Antonin Raﬃn, Ashley Hill, Adam Gleave, Anssi Kanervisto, Maximilian Ernestus,
and Noah Dormann. Stable-baselines3: Reliable reinforcement learning implementa-
tions. Journal of Machine Learning Research, 22(268):1–8, 2021.

[214] Christopher W Moore. The mediation process: Practical strategies for resolving con-

ﬂict. John Wiley & Sons, 2014.

[215] Claude-H´el`ene Mayer.
Springer, 2020.

Intercultural Mediation and Conﬂict Management Training.

[216] Solomon M Hsiang, Marshall Burke, and Edward Miguel. Quantifying the inﬂuence of

climate on human conﬂict. Science, 341(6151):1235367, 2013.

[217] Christopher Cramer. Does inequality cause conﬂict? Journal of International Devel-

opment: The Journal of the Development Studies Association, 15(4):397–412, 2003.

124

[218] Frances Stewart, Douglas Holdstock, and Antonio Jarquin. Root causes of violent
conﬂict in developing countriescommentary: Conﬂict—from causes to prevention? bmj,
324(7333):342–345, 2002.

[219] Gerry O’Reilly. Aligning geopolitics, humanitarian action and geography in times of

conﬂict. Springer, 2019.

[220] Derek J Clark and Kai A Konrad. Asymmetric conﬂict: Weakest link against best

shot. Journal of Conﬂict Resolution, 51(3):457–469, 2007.

[221] J Paul Dunne, Mar´ıa DC Garc´ıa-Alonso, Paul Levine, and Ron P Smith. Managing

asymmetric conﬂict. Oxford Economic Papers, 58(2):183–208, 2006.

[222] Frederick William Lanchester. Aircraft in warfare: The dawn of the fourth arm. Con-

stable limited, 1916.

[223] Dominic DP Johnson and Niall J MacKay. Fight the power: Lanchester’s laws of

combat in human evolution. Evolution and Human Behavior, 36(2):152–163, 2015.

[224] Lucas McDaniel, Erik Talvi, and Brian Hay. Capture the ﬂag as cyber security in-
troduction. In 2016 49th hawaii international conference on system sciences (hicss),
pages 5479–5486. IEEE, 2016.

[225] Fabio Massimo Zennaro and Laszlo Erdodi. Modeling penetration testing with re-
inforcement learning using capture-the-ﬂag challenges: trade-oﬀs between model-free
learning and a priori knowledge. arXiv preprint arXiv:2005.12632, 2020.

[226] Crispin Cowan, Seth Arnold, Steve Beattie, Chris Wright, and John Viega. Defcon cap-
ture the ﬂag: Defending vulnerable code from intense attack. In Proceedings DARPA
Information Survivability Conference and Exposition, volume 1, pages 120–129. IEEE,
2003.

[227] Matthew A Blake, Gerrit A Sorensen, James K Archibald, and Randal W Beard.
Human assisted capture-the-ﬂag in an urban environment.
In IEEE International
Conference on Robotics and Automation, 2004. Proceedings. ICRA’04. 2004, volume 2,
pages 1167–1172. IEEE, 2004.

[228] Michael Novitzky, Paul Robinette, Michael R Benjamin, Danielle K Gleason, Caileigh
Fitzgerald, and Henrik Schmidt. Preliminary interactions of human-robot trust, cog-
nitive load, and robot intelligence levels in a competitive game. In Companion of the
2018 ACM/IEEE international conference on human-robot interaction, pages 203–204,
2018.

[229] Gorka Olalde Mendia, Lander Usategui San Juan, Xabier Perez Bascaran, Asier Bil-
bao Calvo, Alejandro Hern´andez Cordero, Irati Zamalloa Ugarte, Aday Muniz Rosas,
David Mayoral Vilches, Unai Ayucar Carbajo, Laura Alzola Kirschgens, et al. Robotics
ctf (rctf), a playground for robot hacking. arXiv preprint arXiv:1810.02690, 2018.

125

[230] Alexandros Merkouris, Varvara Garneli, and Konstantinos Chorianopoulos. Program-
ming human-robot interactions for teaching robotics within a collaborative learning
open space: Robots playing capture the ﬂag game: Programming human-robot interac-
tions within a collaborative learning open space. In CHI Greece 2021: 1st International
Conference of the ACM Greek SIGCHI Chapter, pages 1–5, 2021.

[231] Eisha Akanksha, Neeraj Sharma, Kamal Gulati, et al. Review on reinforcement learn-
ing, research evolution and scope of application. In 2021 5th international conference
on computing methodologies and communication (ICCMC), pages 1416–1423. IEEE,
2021.

[232] Aur´elien G´eron. Hands-on machine learning with scikit-learn and tensorﬂow: Con-

cepts. Tools, and Techniques to build intelligent systems, 2017.

[233] Christopher JCH Watkins and Peter Dayan. Q-learning. Machine learning, 8:279–292,

1992.

[234] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness,
Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg
Ostrovski, et al. Human-level control through deep reinforcement learning. nature,
518(7540):529–533, 2015.

[235] RICHARD Bellman. Dynamic programming, princeton univ. Press Princeton, New

Jersey, 1957.

[236] Richard S Sutton. Learning to predict by the methods of temporal diﬀerences. Machine

learning, 3:9–44, 1988.

[237] Timothy P Lillicrap, Jonathan J Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez,
Yuval Tassa, David Silver, and Daan Wierstra. Continuous control with deep rein-
forcement learning. arXiv preprint arXiv:1509.02971, 2015.

[238] Kenji Doya. Reinforcement learning in continuous time and space. Neural computation,

12(1):219–245, 2000.

[239] Hado Van Hasselt. Reinforcement learning in continuous state and action spaces.

Reinforcement Learning: State-of-the-Art, pages 207–251, 2012.

[240] Alessandro Lazaric, Marcello Restelli, and Andrea Bonarini. Reinforcement learning in
continuous action spaces through sequential monte carlo methods. Advances in neural
information processing systems, 20, 2007.

[241] Adarsh Sehgal, Hung La, Sushil Louis, and Hai Nguyen. Deep reinforcement learning
using genetic algorithm for parameter optimization. In 2019 Third IEEE International
Conference on Robotic Computing (IRC), pages 596–601. IEEE, 2019.

[242] Felipe Petroski Such, Vashisht Madhavan, Edoardo Conti, Joel Lehman, Kenneth O
Stanley, and Jeﬀ Clune. Deep neuroevolution: Genetic algorithms are a competitive
alternative for training deep neural networks for reinforcement learning. arXiv preprint
arXiv:1712.06567, 2017.

126

[243] Sean Gillies, Casper van der Wel, Joris Van den Bossche, Mike W. Taves, Joshua
Arnott, Brendan C. Ward, et al. Shapely, December 2022. Please cite this software
using these metadata.

[244] Jayesh K Gupta, Maxim Egorov, and Mykel Kochenderfer. Cooperative multi-agent
control using deep reinforcement learning.
In Autonomous Agents and Multiagent
Systems: AAMAS 2017 Workshops, Best Papers, S˜ao Paulo, Brazil, May 8-12, 2017,
Revised Selected Papers 16, pages 66–83. Springer, 2017.

[245] Soummya Kar, Jos´e MF Moura, and H Vincent Poor. Qd-learning: A collaborative
distributed strategy for multi-agent reinforcement learning through consensus + inno-
vations. IEEE Transactions on Signal Processing, 61(7):1848–1862, 2013.

[246] Lucian Bu¸soniu, Robert Babuˇska, and Bart De Schutter. Multi-agent reinforcement
learning: An overview. Innovations in multi-agent systems and applications-1, pages
183–221, 2010.

[247] Thanh Thi Nguyen, Ngoc Duy Nguyen, and Saeid Nahavandi. Deep reinforcement
learning for multiagent systems: A review of challenges, solutions, and applications.
IEEE transactions on cybernetics, 50(9):3826–3839, 2020.

[248] Benjamin Golub and Matthew O Jackson. Naive learning in social networks and the
wisdom of crowds. American Economic Journal: Microeconomics, 2(1):112–149, 2010.

127

APPENDIX A: SUPPLEMENT TO CHAPTER 2

In this supplement we provide a description of our model using the updated Overview,

Design concepts and Details (ODD) protocol [94].

Overview, Design concepts, and Details

Purpose and patterns
The purpose of our model is to produce new, unobserved GPS locations of white-tailed deer
that are consistent with recorded GPS location data of white-tailed deer in the Midwest.
There are three patterns our model aims to reproduce.

1. The ﬁrst pattern is the distribution of jumps (distance between two positions recorded
one after another) should follow a Laplace distribution as observed in our recorded
data.

2. The deer in our data move in one or more distinct spatial locations and occasionally
jump between them. We refer to these spatial locations as basins and the jumps
between them as basin hops. When animals are in one of these spatial locations,
their movements are constrained to a localized region. Therefore the second pattern is
reproducing movements in a single basin.

3. The ﬁnal pattern is reproducing the hops between the basins.

Entities, state variables, and scales
Our model consists of the following three entities.

1. A deer agent that represents a single deer that moves individually in a basin,

2. a herd collective that represents a group of deer that undergo basin hops together, and

3. a world environment containing the herds and deer. The world initializes and runs

each simulation.

The state variables for each of the entities are summarized in Table 8.1, Table 8.2, and Table
8.3.

The spatial scale of our model is a single UTM zone as our GPS data employed the UTM
coordinate system and all data are in a single zone. We generated positions in smaller areas
(approximately 18 km2). We chose this smaller region to examine herd movements that are
close to each other and could overlap. A single time-step in our simulation is 5 hours, since
that was the resolution our GPS was measured.
In our example simulations we ran our
model for 1000 steps (approximately 208 days).

128

Variable name Type and units

x
y

c1x

c1y

bx

by

ﬂoat, dynamic; m
ﬂoat, dynamic; m

ﬂoat, static; unitless

ﬂoat, static; unitless

ﬂoat, static; m−1

ﬂoat, static; m−1

Description
The UTM Easting of a deer.
The UTM Northing of a deer.
Movement model parameter.
Controls size of basin in combination with c1y.
Movement model parameter.
Controls size of basin in combination with c1x.
Laplace distribution scale parameter for UTM
Easting jump distribution.
Laplace distribution scale parameter for UTM
Northing jump distribution.

cov mat

ﬂoat array, static; m2 Covariance matrix of UTM Easting and UTM

Northing data.

Table 8.1: Deer state variables.

Variable name
num deer
cx
cy
num crit jumps data

Type and units
int, static; deer
ﬂoat, dynamic; m
ﬂoat, dynamic; m
int, jumps

xtrds

ytrds

ﬂoat array, static;m

ﬂoat array, static; m

Description
The number of deer in the herd.
UTM Easting of center of basin.
UTM Northing of center of basin.
Scale for Poisson process.
Basin hop trend removed from UTM
Easting data.
Basin hop trend removed from UTM
Northing data.

deer list
jump times
time

Deer list, static; deer object List of deer that are in the herd.
int array, static; time-steps
int, dynamic; time-step

Time-steps which basin hops occur.
Herd’s current time-step.

Table 8.2: Herd state variables.

Variable name Type and units
int, static; herds
int array, static; deer
int tuple list, static; (m,m)
herd list, static; herd object List of herds in simulation.

Description
Number of herds in a simulation.
Number of deer in each herd.
Initial position of each herd.

num herds
herd sizes
herd pos
herd list

Table 8.3: World state variables.

Process overview and scheduling
We designed our model to produce GPS locations of white-tailed deer in the Midwest mea-
sured at approximately 5 hour intervals, as this was the location and time-scale of our data.
There are two movement processes that we wished to capture, localized movement (within

129

basins) and uncommon large movements (basin hops).

Our data does not track every deer in a region, and thus we assumed that groups of
deer (herds) existed and our data tracked one deer in a group. We assumed that all deer in
a shared herd were constrained to the same area (basin), but had independent movements
from each other. When a basin-hop occurred, we assumed all deer in the herd underwent
the large movement together. These processes are captured through the herd’s and deer’s
“move” submodel.

The schedule of each time-step is as follows.

1. The world calls its “update” submodel, which includes:

(a) Each herd calls their “move” submodel, which includes:

i. Each deer calls their “move” submodel, in which they update their individ-
ual positions (x,y) according to the parameters of their herd and a random
component.

ii. The herd checks if a basin-hop should occur on the current time-step. If yes:

A. A basin-hop is sampled and all deer in the herd are moved according to

the sampled jump.

B. The center of the basin is moved according to the sampled jump, updating

cx and cy.

iii. time is incremented by one.

Design concepts

Basic principals

In designing our model we aimed to highlight the importance of EDA in model creation. This
concept is widely applicable to agent-based models of diﬀerent phenomenon. We additionally
wanted to create a model that was able to reproduce three aspects of our data, non-Gaussian
distributed step distributions, localized movement (within basins), and uncommon large
movements (basin-hops). It is common to use Gaussian random walks of some ﬂavor and/or
Levy ﬂights as models for animal movement, but through our EDA we found these models
were not consistent with our data. Therefore, we used our data as a guide to build a model
that was more representative of our data. Namely, our model produces GPS locations of
white-tailed deer using a non-Gaussian random walk in combination with a basin-hopping
model to produce realistic simulated GPS positions. This model could be extended to study
other phenomenon related to deer movement, for example the spread of CWD.

130

Emergence

They key result of our model are the generated GPS locations. Speciﬁcally, jumps (diﬀerences
between adjacent locations) distributions being Laplace distributed, local movements (inside
basins), and uncommon large movements (basin-hops).
In our model, the distribution of
jumps, and two movement patterns are imposed by the deer’s movement update rule. From
the deer’s movement rule, the spatial distribution of deer in a simulation emerge.

Adaptation

Our model has one adaptive behavior, whether or not a herd performs a basin-hop. This
behavior is modeled with indirect objective seeking. There are two rules herds follow to
produce basin-hops. First, herds sample multiple time-steps that they will perform a basin-
hop at. These distributions are trained on real data to produce a number of jumps similar
to the number observed in our data. The second rule is when a herd’s time is equal to one
of the basin-hop times, it samples a basin hop from a kernel density estimate trained on the
observed basin hops, and move all of the deer in its deer list and its center (cx,cy) according
to the sampled jump.

Objectives

Our model only has indirect objective seeking, so there are no objectives.

Learning

The model includes no learning.

Prediction

The adaptive behavior of basin-hopping utilizes implicit prediction. There are many reasons
that deer might disperse, for example ﬁnding a new food source or a dispersal event. Our
model is not attempting to explain how these predictions are actually made, but are modeling
them to produce results that are similar to those observed in our data.

Sensing

Deer are the only agents in this model. They are able to perfectly sense all of their state
variable. Their position (x,y) and movement parameters (c1x, c1y) are used directly to
update their positions, and the remaining variables are used to produce a random component
to their movement. Deer are also able to perfectly sense the center of a the herd they belong

131

to (cx, cy) but not any other herd. They are attracted to the center of the herd they belong
to.

Interaction

Deer do not directly interact with each other, as their inner-basin movements are indepen-
dent. Though deer who belong to a common herd are constrained to the same localized
regions (basins) as they are attracted towards the center of the herd. They also undergo
basin-hops together.

Stochasticity

Our model uses stochasticity in many ways. In the initialization, herds initial center could
be set using a random number generator. Herds are then randomly assigned state variables
by sampling distributions trained on our data to include variability in the herds. Using these
variables, they sample a number of basin-hops they will perform from an Poisson distribution
and when the basin-hops will occur from an exponential distribution. Once the herds have
been created, they populate themselves with deer. The deer’s initial positions are generated
randomly by adding separate random numbers, sampled from a normal distribution with
zero mean and unit variance, to the herds center (cx,cy). These additional random numbers
keep all of the deer from starting at the same location. Deer in diﬀerent herds are randomly
assigned state variables, but deer in a single herd have the same state variables (excluding
x and y)

There is also stochasticity used when the model updates. When a deer moves, it adds a
random component to its new position sampling a distribution constructed from our data.
This keeps deer from making the same moves, and from deer just moving to the center
of the basin they are in. Additionally, when a herd undergoes a basin-hop it samples a
random jump from a distribution constructed from our data. We sample these jumps to add
variability in our model and produce results that diﬀer from our data, but are consistent
with the data.

Collectives

We have one explicit collective in our model, herds. Herds are groups of deer that move in the
same localized areas (basins) and undergo basin hops together. Herds have their own state
variables (see Table 8.2). We included herds as a collective to approximate the movement
of deer that were not tracked in our data. We assumed that the deer that happened to be
tracked had some number of other untracked deer they stayed near. As discussed above the

132

deer in a herd stay in the same basin and undergo basin hops together. It was convenient
to model groups of deer as a collective instead of each individual deer on their own.

Observation

The purpose of our model is to produce new, realistic GPS locations of white-tailed deer.
Therefore the output of our model is the positions of each deer and the herd they belong
to. These data can be used to make density distribution plots of each herd, show the GPS
locations of each herd, and the individual GPS locations of each deer. World has a submodel
“draw” to draw each of the herds GPS locations for all deer. Herds has a submodel “draw”
to draw the GPS locations of each deer in the herd. Finally, deer has a submodel “draw”
that draws its GPS locations.

Initialization
A total number of herds, their sizes, and centers are deﬁned. These parameters along with
distributions of c1x, c1y, bx, by, cov mat, num crit jumps data, xtrds, ytrds are passed to
an instantiation of world.

Each herd is passed their corresponding herd size and center, and all distributions which
were passed to world. The heard then randomly samples the distributions for its state
variables. Then the herd constructs a kernel density estimate from xtrds and ytrds to be
used to sample basin-hops later. Next the herd instantiates the total number of deer passed
from world to deer list. It passes a random initial location (see Sec. 8), and samples the
deers state variables. Each deer in a herd have the same state variables, except x and y.
Next the herd sets the parameter time to zero, and samples the number and times of basin
hops.

For each deer in a herd, it was passes a random initial location and state variables. These

are stored to be used in updating positions.

Input data
The model does not use input data to represent time-varying processes.

Submodels

World

init
simulated.

This submodel instantiates a simulation by spawning all of the herds to be

draw This submodel calls the draw method for all of the herds, which draws the GPS

locations of all the deer.

133

update This submodel calls the “move” method for all of the herds.

Herd

init

This submodel instantiates a herd. It samples the state variables for itself, and
those to be passed to the deer that belongs to the herd. It then spawns all of the deer that
belong to the herd.

get params This submodel samples all of the distributions that are passed to a herd
by world. This is done by generating a uniform random number between zero and the total
number of samples in the distributions. Then the parameters in the position that corresponds
to that random number are returned.

draw This submodel draws the GPS location of all deer in the herd.

make kde This submodel smooths xtrds and ytrds for a herd. Then combines xtrds
and ytrds as a list of two dimensional points. Then each adjacent point is subtracted to
arrive at a list of jumps. Next the average distance between all jumps in list is found. Half
of this average distance was used as the bandwidth for a Gaussian kernel density estimate
(KDE) ﬁt on the list of jump. The resulting KDE can be sampled to generate basin-hops.

move This submodel calls the move function for every deer in the herd. Then it checks
if a basin-hop should occur on the current time-step.
If so it samples a jump from the
KDE produced by “make kde”. This jump is added to the center of the herd (cx,cy) and
every deer’s position in the herd. Finally, time is incremented by one to advance the current
time-step.

lasso times The submodel samples a number of simulated basin-hops from a Poisson
distribution. The argument of the Poisson distribution is the state variable num crit jumps data.
After it generates the number of simulated basin-hops, it samples a time-step for each step
to occur. The cumulative sum of these times is returned.

Deer

init

This submodel instantiates a deer.

draw This submodel plots the GPS location of a simulated deer.

134

rng This submodel uses cov mat to generate normally distributed jumps. It then uses
a copula (this requires the state variables bx and by) to transform the normally distributed
jumps to Laplace distributed ones.

move This submodel updates a deers position. The deer move towards the center of the
basin they belong to, and sample “rng” to add a random component to their new position.

135

APPENDIX B: SUPPLEMENT TO CHAPTER 4

DeGroot Model
In this Supplement, we analyze and adapt a second model, the DeGroot model, [166] to
each of the policies discussed in the main document. The DeGroot model was chosen for
its mathematical transparency, with the goal of assessing its ability to describe the eﬀects
of disinformation-mitigation policies.
In Subsection 8, we review the basic properties of
the DeGroot model.
In the following subsections, the model is then modiﬁed to include
disinformation-mitigation strategies. Section 8, Section 8, and Section 8 contain ﬁgures il-
lustrating the eﬀects of content moderation, education, and counter-campaigns, respectively,
in this model.

DeGroot Model
The model employed in the main document requires moderately intensive computational
eﬀort to explore the wide ranges of anti-disinformation policies, their model implementations
and social graphs.
It would be convenient to have an analytic model of the spread of
disinformation that provides insight and intuition with minimal computational eﬀort. A
model that employs a linear update, contains social-network structure and has an equivalent
of a committed agent is the DeGroot model. [166] We adapt this model to include the three
interventions policies (education, content moderation, and counter campaigns) we consider
in this work and their eﬀects on disinformation dynamics.

In the DeGroot model, we consider a group of N agents connected by a network, which
can be represented by an adjacency matrix A = [aij]. Each agent forms its belief as a
weighted average of the beliefs of its neighbors. Formally, if x(t) = [x1(t), x2(t), . . . , xN (t)]T
represents the beliefs of the N agents at time t, then the belief of agent i at time t + 1 is
given by the discrete-time update

or in matrix notation,

xi(t + 1) =

N
(cid:88)

j=1

tijxj(t),

x(t + k) = Tkx(t).

(8.1)

(8.2)

Here, T = [tij] is referred to as the trust matrix, with tij representing the (constant) weight
that agent i assigns to agent j when updating its belief; in other words, tij quantiﬁes the
extent to which agent i trusts agent j. The trust matrix T is derived from the adjacency
matrix A by normalizing each row to sum to one; i.e., tij = aij/ (cid:80)N
k=1 aik. Starting from

136

some initial set of beliefs x(0), this model repeatedly updates each agent’s belief based on a
weighted sum of all other agents’ beliefs, weighted by levels of trust.

In this framework, we have an interpretable update rule in which tij represents a trust
level that varies from tij = 0 for no trust to tij = 1 for complete trust; the model can
be initialized to these discrete values, or trust levels can be allowed to be continuous. In
Fig. 8.1, we show examples of symmetric, directed graphs with equal weights, and their
corresponding trust matrices. In general, however, the graph represented by A is weighted
and directed: I may listen to you, but you need not listen to me. It is through ﬂexibility in
T that we include committed agents and mitigation strategies in the DeGroot model; in the
context of the DeGroot model, committed agents are referred to as “stubborns”.[248]

Many details of the model dynamics are easily revealed by considering the structure of
T. By construction, T is a right-stochastic matrix and thus has a spectral radius of 1;
this is an important property, as updates are obtained using progressively higher powers of
T, as in (8.2). At least one eigenvalue is 1 and corresponds to a right eigenvector whose
elements are all 1. Eigenvalues associated with left and right eigenvectors are equal for square
matrices; this implies that T has at least one left eigenvector associated with an eigenvalue
of 1. If T can be represented as a directed, weighted graph that is strongly connected (i.e.,
irreducible), then we can apply the Perron-Frobenius theorem that states that there is a
unique left eigenvector v that sums to one and corresponds to an eigenvalue of 1. This
eigenvector determines the ﬁnal consensus belief; speciﬁcally, the ﬁnal state is given by

xﬁnal = vT x(0),

(8.3)

where we see that v represents the long-term inﬂuence of the agents. This equation identiﬁes
v, which is the eigenvector centrality of A, as the relevant parameter for describing the
spread of information, including disinformation. Intuitively, eigenvector centrality measures
a node’s inﬂuence in a network by considering both the number and quality of its connections,
recursively factoring in the inﬂuence of neighboring nodes.
It is also immediately clear
from (8.3) that initial beliefs x(0) contribute to the ﬁnal consensus. A block structure in
T represents independent echo chambers that evolve according to the reducibility class of
each block. It is possible that T is periodic, although some structures (e.g., upper/lower
triangular or symmetric) guarantee that it is aperiodic. We can also examine the second
largest eigenvalue of T to obtain the eigenvalue gap, which reveals the speed of convergence;
larger gaps are associated with faster convergence.

We now extend the DeGroot model to include the three disinformation-mitigation strate-
gies – education, content moderation, and counter campaigns – that we consider in this work

137

Figure 8.1: Example graphs and their corresponding trust matrices. Despite its simplicity,
the DeGroot model contains a complete social network structure as a directed and weighted
graph. Social network graphs are shown in the top row, and their corresponding trust
matrices are shown in the bottom row.

and their eﬀects on the spread of disinformation in the presence of stubborns. We model
these interventions as modiﬁcations to the trust matrix and/or the initial condition and
analyze their impact through the lens of the eigenvalue spectrum and eigenvectors of the
modiﬁed trust matrix. For each mitigation strategy, we begin with a (scale-free) Barabasi-
Albert graph with 100 nodes. Each edge is then converted into two directed arcs that point
in opposite directions, and their weights are determined using social transitivity (described
in the main document). This method of generating directed and weighted edges guaran-
tees that the corresponding trust matrix is irreducible, because the graph is clearly strongly
connected. Opinions have values between 0 and 1; here, we assume that 0 is the truth
and that 1 is disinformation. Values between 0 and 1 indicate a bias towards the truth or
disinformation, with the exception of 0.5, which indicates a completely neutral position.

Stubborns are added to the model by selecting an individual i and setting tii = 1, with
tij = 0 for i (cid:54)= j. Importantly, the addition of stubborns make the trust matrix reducible,
and consensus is not necessarily reached. However, opinions will converge to a value that
depends on the structure of the trust matrix and on agents’ initial opinions. One important
scenario is that in which all stubborns have the same opinion. In such a case, any node that
is reachable from a stubborn will converge to the stubborn’s opinion. Using our method
for generating trust matrices guarantees that the presence of one or more stubborns with
the same opinion, and no stubborns with a diﬀerent opinion, will eventually converge the
network to a consensus of their opinion. Unless otherwise stated, xi(0) = 1 if an individual is
stubborn, and is initially 0 otherwise, in each of our scenarios; i.e., stubborns are committed
to disinformation, and all other agents initially believe the truth. The goal of mitigation
strategies is to minimize xﬁnal when a consensus is reached, or if a consensus is not reached,

138

012345678901234567890123456789then the goal is to push opinions towards the truth or to minimize the eigenvalue gap when
opinions tend towards disinformation.

Education
Skeptical and attentive education strategies are modeled by altering the tij based on the
type of education received by each agent i. In the skeptical education strategy, the values
of the diagonal elements of T of committed agents are decreased; this reﬂects a weaker
commitment to one’s own belief. In the attentive strategy, the tij of neighbors adjacent to
stubborn agents are decreased. To test both education strategies, utilized Barab´asi-Albert
graphs with 100 agents, selecting 10 of these agents to be committed to opinion A, leaving
the remainder uncommitted to opinion B. The edges of these graphs are then transformed
into two directed edges in both directions, weighted using social transitivity as described
in the main document. We generate a trust matrices, {T} by normalizing the adjacency
matrices of these graphs, and then we apply our education mitigation strategies to these
trust matrices.

The skepticism educational strategy reduces diagonal elements tii of the trust matrix
by a factor σ < 1, and the remaining elements corresponding to a committed individual’s
neighbors (determined by the adjacency matrix) are increased by equal parts of σtii, so
that the row still sums to 1. Such a change to T changes the eigenvalue spectrum. Most
importantly, it guarantees that tii < 1, ensuring that the matrix is reducible, and a consensus
opinion is reached at a value less than 1. In Fig. 8.2, we show the consensus opinion versus
σ. For each value of σ, we executed 50 simulations and averaged over initial conditions and
graphs; for each value of σ, the average consensus opinion reached is shown with a black
x, and the value each simulation converged to is shown with a grey dot. When there is
no education (σ = 0), the stubborns’ opinion takes over, and all simulations converge to
disinformation. As the stubborns become more skeptical (σ increases), the average value of
consensus decreases and approaches a constant value of approximately 0.2. The variance in
the consensus value is approximately constant for all values of σ > 0.

The attentive strategy is modeled by increasing elements tij that correspond to non-
stubborn individuals i listening to other non-stubborn individuals j with initial opinions of
zero, by a factor α > 1, and reducing elements ti,k equally such that the rows of T sum to 1.
Such a change biases individuals to listen more to true opinions. The attentive educational
strategy does not aﬀect the opinions of stubborn individuals, which means that eventually,
the population will still converge to the disinformation opinion. However, the attentive
strategy can aﬀect how quickly the model converges, which is measured by the eigenvalue
gap. A smaller eigenvalue gap is associated with a longer convergence time. In Fig. 8.3,
we examine the eigenvalue gap versus α in 50 simulations. The blue line shows the average

139

Figure 8.2: Final consensus opinions versus the education parameter σ. For each value of σ,
gray circles and dark x indicate the value each individual run converged to and the average
ﬁnal opinion over 50 runs, respectively. The ﬁnal opinion decreases toward the truth with
more skepticism education (larger values of σ).

eigenvalue gap in the absence of a mitigation strategy, and the orange line shows the average
eigenvalue gap when long-term education is employed. The 95% conﬁdence intervals for
both are shown as shaded regions. The attentive education mitigation strategy was applied
to the same 50 networks examined in the absence of a mitigation strategy, and thus, both
lines share random ﬂuctuations. When no mitigation is applied, the eigenvalue gap is fairly
constant with respect to α. Attentive education reduced the eigenvalue gap slightly as α
increases, but not substantially.

The results for the skepticism and attentive education strategies we obtained using a
DeGroot model are consistent with what we found when studying such strategies using the
binary agreement model. As we found using that model, reducing people’s commitment to
disinformation can largely counteract the stubborn’s ability to sway individuals’ opinions.
However, solely biasing people towards the truth is not suﬃcient to overcome the inﬂuence
of stubborns.

Content Moderation
In our model of disinformation spread, we ﬁrst create a Barab´asi-Albert graph with 100
agents, selecting 10 of these agents to be committed to opinion A, leaving the remainder
uncommitted to opinion B. The edges of this graph are then transformed into two directed
edges in both directions, weighted using social transitivity as described in the main document.
We create three copies of this graph for distinct analyses. In the ﬁrst copy, we randomly

140

0.000.050.100.150.200.250.20.40.60.81.0opinionindividualaverageFigure 8.3: Eigenvalue gap versus the education parameter α. The blue line shows the
eigenvalue gap when no education is employed, while the orange line shows the eigenvalue
gap for increasing amounts of attentive education. The eigenvalue gap reduces slightly as α
increases meaning the rate that the model converges to the disinformation is slightly slowed.

select n committed nodes to remove, where n varies from 0 to 9. In the second copy, we
systematically remove n committed agents based on their eigenvector centrality. To do this,
we compute eigenvector centrality, remove the node with the highest centrality, recompute
the centrality, and repeat this process until n agents have been removed. The third copy
remains unaltered, serving as a control. For each of these three graphs, we generate a
trust matrix by normalizing the adjacency matrix and then calculate the eigenvalues of the
trust matrices. We also compute the diﬀerence in the two largest eigenvalues, ignoring any
repeated values (e.g., if the eigenvalues are 1, 1, 0.75, then the eigenvalue gap is 0.25).
This entire process is repeated 50 times. For each value of n, ranging from 0 to 9, we
average the 50 runs and plot this average as a solid line, also calculating a 95% conﬁdence
interval to create a band around the solid lines. This methodology allows us to explore the
eﬀects of content moderation through both random and strategic removal of inﬂuential nodes
on the dynamic behavior of the system, as characterized by the eigenvalue gap. Removing
individuals from the graph alters the dominant eigenvector and the eigenvalue spectrum of T,
thereby shifting control away from the removed agents. However, as discussed previously, if
we do not remove all of the stubborns, then the population opinion will eventually converge to
the disinformation opinion. As before, we can then examine the eigenvalue gap to determine
whether the rate of convergence changes.

In Fig. 8.4, we show the eigenvalue gap versus the number of committed individuals

141

0.0000.0250.0500.0750.1000.1250.1500.1750.2000.000.020.040.060.080.100.120.140.16eigenvalue gapno mitigationattentive educationFigure 8.4: Eigenvalue gap versus number of committed agents removed. The eigenvalue
gap versus removing agents randomly is shown in orange, and for removing agents based
on eigenvalue centrality in green. As a baseline we show the eigenvalue gap when not
removing agents in blue. Solid lines are the average of 50 simulations, while shaded regions
are 95% conﬁdence intervals. Removing agents with either strategy reduces the eigenvalue
gap, however, a targeted approach based on eigenvalue centrality reduces the eigenvalue gap
a a faster rate.

removed from the model. Results for removing individuals randomly are shown with a
solid orange line, and results for removing individuals based on their level of inﬂuence are
shown with a solid green line. We also show results for removing no individuals with a solid
blue line, for comparison. The results shown with all three solid lines were obtained using
the same graphs and are averages over 50 simulations. The shaded regions indicate 95%
conﬁdence intervals. When no individuals are removed, the eigenvalue gap is approximately
0.07. When stubborns are removed from the graph with either approach, the eigenvalue gap
tends towards zero. However, when stubborns are removed based on eigenvalue centrality, the
eigenvalue gap converges towards zero faster than when stubborns are removed randomly.
Therefore, in this model, a targeted approach to removing stubborns slows the spread of
disinformation more than does removing individuals randomly.

Counter Campaigns
Counter campaigns are modeled by introducing another group of stubborns who are commit-
ted to the truth. The social networks we consider contain 3 subgroups of agents: 1) agents
spreading disinformation (nA nodes), 2) agents running a counter campaign for truth (nB
nodes), and 3) regular agents comprising the remaining nodes. We model the committed

142

agents as ”stubborn” in the DeGroot model by setting their self-weights to 1. This ensures
that they do not update their opinions over time. The initial opinion vectors for the groups
are xi(0) = 1 if you are committed to A, and 0 otherwise.

In Fig. 8.5 we show many examples of counter-campaigns of various sizes. As in our
previous simulations, we began with a Barab´asi-Albert graph with 100 agents, selecting 10
of these agents to be committed to opinion A, and nB agents to be committed to B, and the
remaining to be uncommitted to B. The edges of this graph are then transformed into two
directed edges in both directions, weighted using social transitivity. Each subplot of Fig. 8.5
shows results for a values of nB between 0 and 10, i.e., for a diﬀerent size of counter-campaign.
In each subplot, each individual’s opinion versus the number of iterations is shown as a grey
line. The committed agents’ opinions can be seen as horizontal lines with values 0 and 1.
Uncommitted individuals’ opinions begin at zero and grow quickly to values between 0 and
1. As discussed previously, when all of the stubborns are committed to disinformation, all
individuals eventually become committed to the disinformation opinion, as can be seen in the
upper left subplot. However, as the size of the counter-campaign grows, uncommitted agents’
opinions converge to opinions that are closer to the truth. When nB = nA, uncommitted
agents’ opinions converge to values near 0.5.

Conclusion
Understanding and combating the spread of disinformation in social networks is a critical
challenge. Through the lens of the DeGroot model, we have examined three potential inter-
ventions and their impacts on the dynamics of belief formation. Our results highlight the
importance of network topology and the distribution of trust in shaping these dynamics. In
comparison with the more complete model in the main document, the DeGroot model allows
for extremely fast exploration of parameter space, mainly related to the cost of diagonalizing
T .

The results we obtained using the DeGroot model resemble the results of the binary
agreement model in the main document. Namely, a skepticism education strategy was by
far the most eﬀective way to combat disinformation, while an attentive one had little to
no aﬀect. In our content moderation policies, we again saw if we do not remove all of the
committed agents, eventually the committed agents will convince everyone of their opinion.
However, in the DeGroot model our targeted approach performed better than our random
approach, however, they both performed similarly. Counter-campaigns seemed to have a
stronger eﬀect in the DeGroot model compared to in the binary agreement model. However
the initial positions of the counter-campaign agents were placed randomly when exploring
the DeGroot model, which possibly gives those agents a further reach than those conisdered
in the binary agreement model.

143

Figure 8.5: Agents opinions versus time in the presence of counter-campaigns. Agent’s
opinions versus the number of iterations of the DeGroot model are shown in grey for various
sizes of counter-campaigns. For each size of counter-campaign (nB), the size of the minority
committed to the disinformation is held constant (nA = 10). As the counter-campaign grows,
uncommitted agents converge to opinions closer to the truth. When the counter-campaign
is equal in size to the minority committed to disinformation, uncommitted agents converge
to opinions near .5.

144

The DeGroot model allowed us to quickly explore some of our mitigation strategies, but
its simplicity limits the applicability of results derived from using it. Utilizing the binary
agreement model in the main document allowed us to incorporate more complex attributes,
such as nonlinearities. Nonetheless, further research is needed to develop more sophisticated
models and strategies for combating disinformation in social networks.

Content Moderation Figures
This section contains all of the results for content moderation applied to various artiﬁcial
social networks. Each plot corresponds to the graph type listed in the caption. Results for
removing nodes based on their level of inﬂuence are shown with solid lines, and those for
removing nodes randomly are shown with dashed lines. The colors of the lines represent the
percentage of the total nodes that were removed. Each opaque line shows the average of 5
simulations, and each transparent line shows results for a single simulation.

Figure 8.6: Barbell graph.

145

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removedFigure 8.7: Complete graph.

Figure 8.8: Erdos-Renyi graph with p = .02.

146

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removed0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removedFigure 8.9: Erdos-Renyi graph with p = .12.

Figure 8.10: Erdos-Renyi graph with p = .25.

147

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removed0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removedFigure 8.11: Grid graph.

Figure 8.12: Watts-Strogatz random graph with k = 8 and p = 1.

148

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removed0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removedFigure 8.13: Watts-Strogatz random graph with p = 1 and k = 48.

Figure 8.14: Watts-Strogatz small world graph with p = 1 and k = 100.

149

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removed0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removedFigure 8.15: Watts-Strogatz small world graph with p = .5 and k = 8.

Figure 8.16: Watts-Strogatz small world graph with p = .5 and k = 48.

150

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removed0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removedFigure 8.17: Watts-Strogatz small world graph with p = .5 and k = 100.

Figure 8.18: Barabasi-Albert graph with m = 4.

151

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removed0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removedFigure 8.19: Barabasi-Albert graph with m = 4.

Figure 8.20: Barabasi-Albert graph with m = 24.

152

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removed0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removedFigure 8.21: Barabasi-Albert graph with m = 50.

153

0.030.040.050.060.070.080.090.100.110.120.13pa0.00.20.40.60.81.0nBnone removed0.25% removed0.5% removed1.0% removed2.0% removed2.5% removedEducation Figures
This section contains all of the results for applying early and late education to various
artiﬁcial social networks. Each plot corresponds to the graph type listed in the caption. We
show the eﬀects of early education with an orange dashed line, late education with a green
dotted line, a combination of both with a solid blue line, and no education with a black
dot-dashed line. Each opaque line shows the average of 30 runs, and transparent lines show
results for individual runs.

Figure 8.22: Barbell graph.

154

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattentionFigure 8.23: Complete graph.

Figure 8.24: Erdos-Renyi graph with p = .02.

155

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattention0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothearlylateFigure 8.25: Erdos-Renyi graph with p = .12.

Figure 8.26: Erdos-Renyi graph with p = .25.

156

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothearlylate0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattentionFigure 8.27: Grid graph.

Figure 8.28: Watts-Strogatz random graph with p = 1 and k = 8.

157

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattention0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothearlylateFigure 8.29: Watts-Strogatz random graph with p = 1 and k = 48.

Figure 8.30: WSR25.

158

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothearlylate0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattentionFigure 8.31: Watts-Strogatz small world graph with p = .5 and k = 8.

Figure 8.32: Watts-Strogatz small world graph with p = .5 and k = 48.

159

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothearlylate0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothearlylateFigure 8.33: Watts-Strogatz small world graph with p = .5 and k = 100.

Figure 8.34: Barabasi-Albert graph with m = 4.

160

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattention0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattentionFigure 8.35: Barabasi-Albert graph with m = 24.

Figure 8.36: Barabasi-Albert graph with m = 50.

161

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattention0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBnobothskepticismattentionCounter-Campaign Figures
This section contains all of the results for applying counter-campaigns to various artiﬁcial
social networks. Each plot corresponds to the graph type listed in the caption. As before,
we show the results for no intervention with a black dot-dashed line. We show the eﬀects of
a small counter-campaign (pb = .05) with a solid blue line and of a large counter-campaign
(pb = .15) with a dashed orange line. Each opaque line shows the average of 10 simulations,
and transparent lines show the results for individual runs.

Figure 8.37: Barbell graph.

162

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.15Figure 8.38: Complete graph.

Figure 8.39: Erdos-Renyi graph with p = .02.

163

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.150.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.15Figure 8.40: Erdos-Renyi graph with p = .12.

Figure 8.41: Erdos-Renyi graph with p = .25.

164

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.150.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.15Figure 8.42: Grid graph.

Figure 8.43: Watts-Strogatz random graph with p = 1 and k = 8.

165

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.150.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.15Figure 8.44: Watts-Strogatz random graph with p = 1 and k = 48.

Figure 8.45: WSR25.

166

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.150.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.15Figure 8.46: Watts-Strogatz small world graph with p = .5 and k = 8.

Figure 8.47: Watts-Strogatz small world graph with p = .5 and k = 48.

167

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.150.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.15Figure 8.48: Watts-Strogatz small world graph with p = .5 and k = 100.

Figure 8.49: Barabasi-Albert graph with m = 4.

168

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.150.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.15Figure 8.50: Barabasi-Albert graph with m = 24.

Figure 8.51: Barabasi-Albert graph with m = 50.

169

0.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.150.000.050.100.150.200.250.300.35pa0.00.20.40.60.81.0nBpb=0pb=.05pb=.15