# Modelling host population support for combat adversaries

Mathew Zuparic a, Sergiy Shelyagb, Maia Angelovab, Ye Zhub and Alexander
Kalloniatis a a Defence Science and Technology Group, Canberra, ACT 2600
Australia;
b Data to Intelligence Research Centre, Deakin University, Geelong, VIC 3220
Australia.

###### Abstract

We consider a model of adversarial dynamics consisting of three populations,
labelled Blue, Green and Red, which evolve under a system of first order
nonlinear differential equations. Red and Blue populations are adversaries and
interact via a set of Lanchester combat laws. Green is divided into three sub-
populations: Red supporters, Blue supporters and Neutral. Green support for
Red and Blue leads to more combat effectiveness for either side. From Green’s
perspective, if either Red or Blue exceed a size according to the capacity of
the local population to facilitate or tolerate, then support for that side
diminishes; the corresponding Green population reverts to the neutral sub-
population, who do not contribute to combat effectiveness of either side. The
mechanism for supporters deciding if either Blue or Red are too big is given
by a logistic-type interaction term. The intent of the model is to examine the
role of influence in complex adversarial situations typical in counter-
insurgency, where victory requires a genuine balance between maintaining
combat effectiveness and support from a third party whose backing is not
always assured.

††articletype: ARTICLE TEMPLATE

Keywords Lanchester model; Volterra-Lotka model; influence modelling

##  1 Introduction

The intent of this work is to explore the mutual influence between non-
combatants and combatants in warfare. To enable this we devise a scenario
shown in Figure 1 which contains three distinct populations, labeled Blue, Red
and Green. The Blue and Red forces are homogeneous and in direct opposition
with each other. Green represents a non-homogeneous civilian population in
proximity to the Blue and Red forces, divided into three sub-populations. Sub-
populations $\gamma$ and $g$ support the Blue and Red forces, respectively.
The third sub-population $\Gamma$ remains completely neutral, offering no
support for either Blue or Red. The heterogeneity of Green is similar to the
model proposed in (Atkinson et al. 2012), where the authors divided their Blue
and Red populations into regions of supporters and opponent sympathisers.
Arguably, this is a situation often encountered in counter-insurgency where
government forces (in recent history, supported by external actors) contest
opponents in a space where members of a population between the two may be
swayed in one direction or the other.

Figure 1: Diagram showing Blue, Red and Green populations, and their
interactions.

The model dynamics reflects the arrows which bridge the boundaries of the
populations in Figure 1. Blue and Red affect each other negatively via a
Lanchester-type interaction. The support offered to Red and Blue from the sub-
populations $g$ and $\gamma$ respectively, is determined by the magnitude of
the sub-populations. Additionally, the model allows for dynamic movement
amongst the three Green sub-populations. The movement from one sub-population
to another in this model represents a shift in the Green population’s cultural
perception of the conflict between Blue and Red, akin to the recent
application of the Axelrod cultural model by (González-Avella et al. 2014).
The direction and rate of the flow between the Green sub-populations is
determined by a logistic-type expression; if either Blue or Red become too big
the support they derive from their corresponding Green sub-populations begins
to wane.

###  1.1 Mathematical modelling of adversarial settings

An important milestone of adversarial mathematical modelling can be traced
back to Verhulst’s (1838) formulation of the logistic equation:

$\dot{P}=\beta P\left(1-\frac{P}{k}\right),\;\;P(0)=P_{0},$ (1)

for population $P\in\mathbb{R}$. In Eq.(1), the tension between the linear
growth term with growth rate $\beta\in\mathbb{R}_{+}$, and the quadratic decay
term with carrying capacity $k\in\mathbb{R}_{+}$, ensures population
equilibrium at the carrying capacity value. The tension between Eq.(1)’s
linear growth and nonlinear decay remains a central theoretical mechanism for
the modelling of single species (Richards 1959), and multi-species dynamics
(Royama 1971).

The force attrition model of Osipov (1995) and Lanchester (1916) is also a
popular mathematical model to understand combat outcomes. The model has two
typical variants given by

$\displaystyle\begin{split}&\dot{B}=-\beta_{R}R,\;\;\dot{R}=-\beta_{B}B,\;\;\textrm{square
law/direct fire},\\\
\textrm{or}\;\;&\dot{B}=-\beta_{R}BR,\;\;\dot{R}=-\beta_{B}BR,\;\;\textrm{linear
law/indirect fire},\end{split}$ (2)

where $B$ and $R$ correspond to the magnitudes of the opposing Blue and Red
forces with initial conditions $\\{B(0),R(0)\\}=\\{B_{0},R_{0}\\}$, and the
parameters $\beta_{B}$ and $\beta_{R}$ represent their respective rates of
effectiveness. Eq.(2) was devised as a simple model to understand the
interplay between total force numbers and fighting effectiveness in static
environments (MacKay 2006). Generalisations of the Lanchester equations
include Deitchman’s guerrilla warfare model (1962), which mixes the linear and
square Lanchester laws to account for the effectiveness of an opponent
employing unconventional tactics, and MacKay’s mixed forces model (2009; 2012)
for a general number of $N$ Blue forces vs $M$ Red forces.

Previous studies have developed the Lanchester model for more realistic
lethality rates (Taylor and Parry 1975; Taylor and Brown 1976); applied
optimisation and game theoretic considerations (Taylor 1977, 1979; Lin et al.
2014; MacKay 2015); and compared model outputs to the real world data sets of
World War II (Lucas and Turkes 2004; Hung et al. 2005), the Vietnam War
(Schaffer 1968), and the multi-adversary complexity of the Syrian civil war
(Kress et al. 2018). More recently, (Kalloniatis et al. 2021) proposed a fully
networked Lanchester model and established the emergence of behaviour
reminiscent in manoeuvre warfare theory. Building on the theme of adding
heterogeneity to Eq.(2), recently (Keane 2011) and (González and Villena 2011)
built on the work of (Protopopescu et al. 1989) and added a spatial reaction-
diffusion component to Lanchester dynamics, which enabled concepts such as
movement, terrain and perception to be reflected in the model.

The predator-prey model of Lotka (1925) and Volterra (1928) is another
adversarial model this work takes inspiration from. The simplest two-species
instance of the model is given by,

$\dot{N}=\beta^{(N)}_{g}N-\beta^{(N)}_{d}NP,\;\;\dot{P}=\beta^{(P)}_{g}NP-\beta^{(P)}_{d}P,$
(3)

describing the interaction between a prey-species $N$, and predator-species
$P$. Since its original inception, the highly idealised system in Eq.(3) has
been generalised to include more realistic predation (Holling 1966), a general
number of interacting species (Hening and Nguyen 2018; Hening et al. 2018),
and stochastic terms to simulate environmental effects (Arató 2003).

###  1.2 Exploring the role of influence in adversarial settings

Of particular interest to this work is the role that influence plays in
adversarial settings. This work is inspired by past OR studies which have
applied ODEs to model insurgencies and cooperative dynamics. For example,
(MacKay 2015) combined the Lanchester, Deitchman and Richardson models in an
attempt to demonstrate that state powers cannot deal with insurgencies solely
by force; (Kress and Szechtman 2009; Kress and MacKay 2014) generalised the
Deitchman model to explore the role that intelligence plays in counter-
insurgency operations; (Syms and Solymar 2015) combined the Lanchester and
logistic equations to model regime changes arising from insurgencies; and
(McLennan-Smith et al. 2019, 2020) generalised Lotka-Volterra models to
include symbiotic interactions to model adversarial dynamics in the vicinity
of a homogeneous non-combatant humanitarian agencies.

This current work is situated within the framework of multi-party Lanchester
type models, such as (MacKay 2009, 2012; Kress et al. 2018; Kalloniatis et al.
2021). As alluded, this approach is particularly relevant to counter-
insurgency modelling, though not wholly bound by that context. We construct
two relatively simple models which offer a means to explore the mutual
influence between combatants and non-combatants in warfare. This is achieved
by building upon the Lanchester-inspired ODE approaches cited in this work,
and beyond, by adding inhomogeneity via logistic-Lotka-Volterra interactions
to the non-combatant population. In particular, we exploit the ecology lens
which has been richly discussed for counter-insurgency for many years, though
typically in qualitative analogical terms (Drapeau 2008; Johnson 2008;
Zimmerman 2016; Kilcullen 2020).

In contrast to (McLennan-Smith et al. 2019, 2020) which also exploits
ecological mechanisms to model complex conflict environments, the present work
presents the context where the civilian population may influence either force
in their combatant role, rather than as humanitarian contribution. Naturally,
the modern battlefield has both elements, as evidenced in concepts such as the
three-block war (Dorn and Varey 2009).

###  1.3 Outline of the paper

In section 2 we present two instances of the model, one which conserves the
total Green population over time — labelled the supporter model — and one
which lacks the conservation of the Green population — labelled the
contributor model. To each model we offer example trajectories which inform
our expectations for model behaviours. In section 3 we perform parameter
sweeps which expose a number of interesting behaviours demonstrated by each of
the models and discuss the likely reasons behind these outcomes. We also
perform semi-analytic approximations of the model in a number of regimes in
order to understand mechanisms behind these behaviours. Finally in section 4
we discuss various implications stemming from outcomes and behaviours
witnessed from the model and flag the potential for further work.

##  2 The models

The two models we present in the following may be seen as a segregation of a
more complex unified model, which we separate in order to understand the
dynamics peculiar to each. The overall framework is that using an ecology lens
to understand counter-insurgency, such as Figure 2 of (Drapeau 2008). Here, a
spectrum of discrete states between ‘government’ forces, here seen as Blue,
and ‘insurgent’ combatants, or Red is decomposed: government infrastructure
that provides resource into Blue, government supporters, a broad political
space occupied by the population, then underground supporters of the
insurgency, and an insurgent auxiliary that can be drawn into the conflict as
Red combatants. Thus we consider two segmentations of this spectrum: one
between supporters and combatants, and the other between contributors and
combatants. These each require quite different types of non-linearities with
implications we seek to understand.

###  2.1 Model I: Green supporters

Consider three populations, Blue, Green and Red, labeled $B,G$ and $R$
respectively. $B$ and $R$ are adversaries, detrimentally affecting each other
akin to the directed-fire Lanchester-type interaction. The population $G$, is
composed of three sub-populations, $G=g+\gamma+\Gamma$. Sub-populations $g$
and $\gamma$ are Red and Blue supporters, respectively. Sub-population
$\Gamma$ maintains neutrality and support neither side.

The Lanchester interaction between Blue and Red is given by,

$\displaystyle\begin{split}\dot{B}=-k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}f(g,\gamma)R\Theta(B),\;\;\dot{R}=-k^{\textrm{\tiny{L}}}_{\textrm{\tiny{B}}}f(\gamma,g)B\Theta(R),\\\
\textrm{where}\;\;f(x,y)=\frac{x}{y+1},\end{split}$ (4)

with the lethality coefficients
$\\{k^{\textrm{\tiny{L}}}_{\textrm{\tiny{B}}},k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}\\}\in\mathbb{R}_{+}$.
Both expressions in Eq. (4) are modulated by the dimensionless term $f$, which
factors the impact of both supporter populations, $g$ and $\gamma$. The first
argument in $f$ is the size of the supporting population, and serves to
increase the supported force’s combat effectiveness against their adversary.
For example, this may represent moral support which enables application of
maximal combat power, or the utility of the supporting population’s
infrastructure. Correspondingly, the second argument in $f$ is the size of the
detracting population, decreasing the aforementioned force’s effectiveness
against the adversary. The term ‘1’ in the denominator of $f$ represents a
‘standing population’ that the detractors wish to maintain. If replaced by a
number greater or less than than unity, the detractor population finds it
easier or harder to influence the outcome of combat, respectively. Obviously,
a choice close to zero would potentially lead to pathological behaviour in the
model due to the appearance of unrealistically large values. Furthermore, the
quantity $\Theta$ in both expressions in Eq. (4) is a sufficiently smooth
numerically stable Heaviside-step-like function which ensures that the
populations of Blue and Red do not become negative. This function requires two
parameters: the steepness of the transition and the offset from zero at which
an entity is deemed extinct. For the remainder of this work we use

$\Theta(x)=\frac{1}{2}\left\\{\tanh\left[\alpha\left(x-\frac{4}{\alpha}\right)\right]+1\right\\},$
(5)

where $\alpha=10^{6}$ sets the scale for both the steepness and the offset,
chosen sufficiently large for a sharp transition, and the offset, up to this
scale, is set at ‘4’. Intuitively, for sharp transitions the resultant
dynamics will not be acutely sensitive to small variations in this offset.

The dynamics amongst the sub-populations of Green are determined by,

$\displaystyle\begin{split}\dot{g}=g\Gamma
R\left(1-\frac{R}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}\Theta(g)\right),\\\
\dot{\gamma}=\gamma\Gamma
B\left(1-\frac{B}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}}\Theta(\gamma)\right),\\\
\dot{\Gamma}=-g\Gamma
R\left(1-\frac{R}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}\Theta(g)\right)-\gamma\Gamma
B\left(1-\frac{B}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}}\Theta(\gamma)\right),\end{split}$
(6)

where the inclusion of the $\Theta$ Heaviside terms ensure that the Green sub-
populations remain positive. We remark that there is conservation of Green’s
total population i.e. $\dot{G}\equiv\dot{g}+\dot{\gamma}+\dot{\Gamma}=0$.

The parameters
$\\{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}},k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\\}\in\mathbb{R}_{+}$
in Eq. (6) are carrying capacities associated with the corresponding logistic
terms. The logistic terms are the mechanisms employed to increase/decrease the
sub-population sizes. Here the carrying capacity plays a significant role,
well recognised in insurgency modelling (Drapeau 2008). However, the
limitations of one side or the other to ‘carry’ support go beyond purely
financial (Syms and Solymar 2015) or terrain (Johnson 2008) factors, nor even
humanitarian (Williamson 2011) but also socio-cultural (Metz and Millen 2004).
This includes the social and tactical discipline of the respective combatants.
Thus, if the supported force ($R$ or $B$) becomes larger than the carrying
capacity ($k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$ or
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$) of their respective supporting
Green sub-population ($g$ or $\gamma$), then the rate of change of the sub-
population will become negative. This signals a rescinding of support,
interpreted as non-combatants responding to negative consequences of military
forces engaging in close proximity. Such consequences include collateral
damage, or cultural misunderstanding. Thus a higher carrying capacity value
($k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$ and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$) corresponds to a more tolerant
non-combatant population, and/or a military force which prioritises the best
interests of the non-combatant population. The classic example of this in 21st
century experience is the Anbar Awakening which presaged the 2007 surge of US
forces in Iraq (Zimmerman 2016). These forces would possess adequate training,
with doctrine which encourages military leaders to engage with community
leaders and facilitate their role in decision-making (Metz and Millen 2004;
Petraeus 2009). Evidently, different contexts will have significantly
different population carrying capacities for a Blue force, contrasting the US
occupations, respectively, of Germany after World War II, and Iraq after the
Second Gulf War.

In (Kress and Szechtman 2009) the authors proposed a Lanchester-insurgency
model which considered the effect of a homogeneous and static Green
population. The inclusion of the logistic terms in Eq.(6) enable a
generalisation of the model in (Kress and Szechtman 2009) by considering the
effects of a non-homogeneous Green which can both offer and rescind support to
both forces.

Figure 2: Examples of solutions for Model I. The red and blue curves
correspond to $R$ and $B$ populations, respectively. The dotted green curve
corresponds to $\gamma$ (support for $B$), the solid curve - to $g$ (support
for $R$), the dashed green curve - to the neutral population $\Gamma$, and the
dash-dotted green curve to the total $G$. The initial conditions are $B(0)=2$,
$R(0)=1$, $g(0)=1$, $\gamma(0)=1$, $\Gamma(0)=3$ and $B(0)=1.5$, $R(0)=1.5$,
$g(0)=1$, $\gamma(0)=2$, $\Gamma(0)=3$, for the top and bottom rows,
respectively. The lethality coefficients are always
$k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}=1$ and
$k^{\textrm{\tiny{L}}}_{\textrm{\tiny{B}}}=1$. The carrying capacities are:
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=2$,
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=2$ (left column) and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=1$,
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1$ (right column).

Examples of solutions for Model I are shown in Figure 2. For all cases,
demonstrated in the figure, despite the lethality coefficients being the same
for Blue and Red populations, the outcome is not related to the initial
populations, but strongly depends on the carrying capacity values. Namely, for
the initially different Blue and Red populations (Blue is greater than Red;
top row in Figure 2), Red wins for
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=1$ and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1$, while for larger capacities
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=2$ and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=2$, Blue wins. This is due to the
support offered by Green (solid green curve).

###  2.2 Model II: Green contributors

In the second model instance, the sub-populations $g$ and $\gamma$ play the
role of combat contributors to Red and Blue, respectively. As with the
previous supporter model, sub-population $\Gamma$ maintains neutrality. The
Lanchester-type interaction between Blue and Red is given by,

$\displaystyle\begin{split}\dot{B}=-R\Theta(B)+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\gamma\Theta(\gamma),\;\;\dot{R}=-B\Theta(R)+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}g\Theta(g).\end{split}$
(7)

Eq. (7) differs from the interactions presented in Eq. (4) by placing the sub-
populations of $\gamma$ and $g$ on the same level as the adversaries $B$ and
$R$, rather than as arguments into the modulating the dimensionless factor
$f$. Sub-populations $\gamma$ and $g$ now act as sources, and are employed to
solely decrease the effectiveness of forces $B$ and $R$, respectively. Thus,
the meaning of the parameters
$\\{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}},k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}\\}$
has changed from the previous case, being referred to Blue and Red’s transfer
coefficients, respectively. Additionally, compared to Eq.(4), we have set the
lethality coefficients to unity in this instance, which can be generalised in
the future. Again, the inclusion of the multiplicative $\Theta$ terms ensures
the positivity of population values. In order to reflect the different way
that $\gamma$ and $g$ interact in the adversarial setting in Eq. (7), we
change the dynamics amongst the three sub-populations of Green to be non-
conservative via,

$\displaystyle\begin{split}\dot{g}=g\Gamma
R\left(1-\frac{R}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}\Theta(g)\right)-k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}g\Theta(g),\\\
\dot{\gamma}=\gamma\Gamma
B\left(1-\frac{B}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}}\Theta(\gamma)\right)-k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\gamma\Theta(\gamma),\\\
\dot{\Gamma}=-g\Gamma
R\left(1-\frac{R}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}\Theta(g)\right)-\gamma\Gamma
B\left(1-\frac{B}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}}\Theta(\gamma)\right).\end{split}$
(8)

As with Eq. (6) the constants
$\\{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}},k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\\}\in\mathbb{R}_{+}$
play the role of carrying capacities, with their function and implied meaning
unchanged from the first model discussed in section 2.1. Furthermore, the
defining equations for $g$ and $\gamma$ now contain additional sink terms,
balancing the source terms in the defining equations for Blue and Red. Thus,
the total rate of change of $G$ is now given by
$\dot{G}=-k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\gamma\Theta(\gamma)-k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}g\Theta(g)$.

Figure 3: Example outputs of Model II for various parameter choices. Colour
meanings for trajectories are same as in figure 2. The initial conditions are
$B(0)=20$, $R(0)=20$, $g(0)=10$, $\gamma(0)=10$, $\Gamma(0)=10$ and same for
all given examples. The transfer coefficients are always
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}=1$ and
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}=1$. The carrying capacities are:
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=10$,
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=10$ (top-left plot);
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=1$,
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1$ (top-right plot);
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=1$,
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=10$ (bottom-left plot);
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=10$,
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=100$ (bottom-right plot).

This model instance focuses on the same initial conditions for populations $B$
and $R$, and $\gamma$ and $g$, in order to understand the effect that Green
has on the outcome for near-peer adversaries. Additionally, varying initial
conditions in the contributor model matters less as $B$ and $R$ can now grow
in population, resulting in both populations potentially sampling more of the
phase space for each trajectory. Future work concerning this model may address
different initial conditions more thoroughly. Examples of solutions for this
model instance are shown in Figure 3. Larger initial conditions were chosen to
more easily demonstrate details in trajectories. For all panels the transfer
coefficients are the same with
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}=k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}=1$.
Outcomes and trajectories differ greatly however, strongly depending on the
values of the carrying capacities. Focusing on the top row for equal carrying
capacities the result is a draw; though this stalemate is achieved through an
oscillatory stage, whose temporal extent depends on the carrying capacities.
For large
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=10$
(top-left), the oscillatory stage is short, and the Green population decreases
faster than the Blue or Red populations. For small
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1$
(top-right), the Red and Blue populations quickly decrease to almost zero,
while the Green population keeps its oscillatory decrease until exhausted. We
will address the nature of oscillations in the solution in Section 3.

For different carrying capacities, the outcome depends on whether one of the
carrying capacities is greater than the corresponding Red or Blue population.
For $k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=1$ and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=10$ (bottom-left; both carrying
capacities are less than $B(0)$ and $R(0)$), the Red and Blue populations
initially decrease until Blue reaches the threshold level. Then, Blue
population receives support from the Green population, which consequently
leads to Blue’s victory. However, for
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=10$ and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=100$ (bottom-right;
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}>B(0)$), Blue gets immediate support
from the Green population, which leads to Blue victory. The behaviour of the
Green population is very different in the two latter cases. For small carrying
capacities, the Green population remains constant before and after the
critical phase at $t=1$, and a small decrease in Green population leads to the
Blue victory. For large carrying capacities, Green contributes to Blue over
the extent of the solution and keeps contributing to Blue even after Red
population reaches zero at $t=0.7$.

##  3 Numerical analysis and semi-analytic approximations

In order to understand model behaviours for both model iterations, we offer
more substantial numerical analysis which highlight differences due to various
parameters changes. Furthermore, in key instances we solve appropriate model
approximations which enable a more in-depth understanding of interesting model
behaviours.

###  3.1 Supporting case: parameter sweeps

For the Green-supporting scenario given by Eqs.(4,6), we offer phase plots in
Figure 4 which detail final values (taken at $t=50$) for the populations over
a range of $k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}\in(0,20)$ (Red’s
lethality coefficient) and $k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\in(0,2)$
(carrying capacity dictating $g$’s capacity to support Red). Thus, in Figure 4
we are examining model outcomes of Red changing its parameter inputs — combat
effectiveness, and it’s ability to maintain support from $g$. In all four
panels the initial conditions are $B_{0}=R_{0}=1.5$, $g_{0}=1$, $\gamma_{0}=2$
and $\Gamma_{0}=3$.

Figure 4: Parameter scans for Model I, $B(t=0)=1.5$, $R(t=0)=1.5$, $g(t=0)=1$,
$\gamma(t=0)=2$, $\Gamma(t=0)=3$,
$k^{\textrm{\tiny{L}}}_{\textrm{\tiny{B}}}=1$,
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1$. Colours indicate the final
population values (at $t=50$) over a range of
$k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}$ and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$, with blue and red signifying $B$
and $R$ victories in the top left panel. The white contour in the top-left
panel separates the R and B victory regions.

The top left panel of Figure 4 details the final value of the expression $R-B$
for Eq.(4). Positive values (Red wins) are indicated by shades of red, and
negative values (Blue wins) indicated by shades of blue. White contours
separate the regions between Red and Blue victories. In the region
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}>1$ Red always wins engagements.
Correspondingly, for low values of
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}<0.2$ Blue is always the victor. In
the region $0.2<k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}<1$ the outcome is
more complicated due to the presence of two white contours (signifying a
stalemate) separating the regions of Red and Blue winning. Thus in this region
there is a thin slice between the majority of Blue victories where Red emerges
victorious.

Contrast to this, on the bottom-left panel for Red contributors $g$, the
regions where $g$ are non-zero correspond to the regions where Red wins the
engagement in the top left diagram. This includes the appearance of the thin
slice between the Blue victories in the region
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\in(0.2,1)$. This slice also exists
as non-zero $g$ values in the bottom-left panel. Nevertheless, there is a
notable half-oval shape in the parameter values
$k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}\times
k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\in(10,20)\times(1,1.5)$, absent from
the top-left diagram, indicating a region of zero $g$-values, even though Red
wins the engagement in this region. Trajectories in this region reveal that
Blue reaches zero population before Red descends below
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$-value — i.e. Red does not have the
ability to sustain support from Green. Thus in this region $g$’s population
declines after the engagement ends. Just below this region in the
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$ direction Red still obtains
victory, but the final values of $g$ are non-zero. Trajectories in this region
show that although $g$-values initially begin decreasing, once Blue loses the
engagement Red population is below
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$-value, meaning that $g$ increases
after the engagement has ended. Decreasing
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$ values there is an enclosed region
where Blue obtains victory, in addition to small (or zero) $g$-values at the
final time. Trajectories in this region reveal that although Blue starts the
engagement against Red poorly, because
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}<k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$,
and hence Red’s ability to sustain Green support is less than Blue’s, Blue
eventually wins the engagement due to experiencing much more support.

Focusing on lower $k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$ values, we reach
the aforementioned strip of red surrounded by Blue victories. In this region
Red wins the engagement in combination with large $g$-values at the final
time. Trajectories in this region reveal that Blue is actually ahead in
engagements for the majority of the time. However, as Red’s population value
for the majority of time is less than
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$, whereas Blue’s is greater than
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$, Red is able to maintain support
from $g$ for a significant amount of time, whereas Blue’s support from
$\gamma$ dwindles away. Red’s growing support from Green means that Blue’s
combat effectiveness eventually plummets, leading to Red victory. In the
region for lowest $k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$ values which
show Blue victory in the top panel, there are zero $g$ and $\gamma$ values. In
this region, both Red and Blue completely lose support from their respective
Green populations. Nevertheless, due to Blue’s ability to sustain Green
support being greater than Red’s, Red has less (though non-zero) force than
Blue by the final time.

###  3.2 Supporting case: analytic considerations

Using approximations it is possible to gain an analytical understanding of the
mechanisms causing the bottom white contour in the top-left plot of Figure 4,
below which Green has completely rescinded support from both Red and Blue. To
enable this we assign
$k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{L}}}_{\textrm{\tiny{B}}}$
and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$
in Eq.(6), and consider the case $\gamma=g$, thus focusing on the case where
Green supports both forces equally. Furthermore, if we assume that the
starting populations of Blue and Red are equal (i.e $B_{0}=R_{0}$), then this
forces $B=R$ in Eq.(4). The significance of this assumption is that we focus
on what may be deemed the worst-case scenario of stalemate where conflict is
drawn out with devastating and exhausting impacts on the people and resources
of both sides of the conflict.

Thus the system under consideration becomes

$\dot{B}=-k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}\frac{g}{1+g}B,\;\;\dot{g}=g(G_{0}-2g)B\left(1-\frac{B}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}\right),$
(9)

where the initial condition for $g$ is bounded by $0<g_{0}<G_{0}/2$. Eq.(9)
reveals two important facts:

  * •

the right hand side of the equation for $\dot{B}$ is always less than or equal
to zero

  * •

if $B_{0}\leq k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$, then the right hand
side of the equation for $g$ will remain positive, ensuring that $B\rightarrow
0$ after sufficient time.

We solve Eq.(9) to obtain the relationship between dynamic variables $B$ and
$g$ and determine the remaining conditions for
$B_{0}>k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$ which enable $B\rightarrow
0$ after sufficient time. Dividing $\dot{B}$ by $\dot{g}$ in Eq.(9) reveals

$\displaystyle\begin{split}&\frac{dB}{dg}=-\frac{k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}}{(1+g)(G_{0}-2g)\left(1-\frac{B}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}\right)},\\\
\Rightarrow\;\;&(B-B_{0})(B+B_{0}-2k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}})\\\
&=\frac{2k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}{G_{0}+2}\ln\left[\left(\frac{1+g}{G_{0}-2g}\right)\left(\frac{G_{0}-2g_{0}}{1+g_{0}}\right)\right].\end{split}$
(10)

Expressing $g$ as a function of $B$, Eq.(10) becomes

$\displaystyle\begin{split}&g=\frac{G_{0}h(B)-1}{1+2h(B)},\\\
\textrm{where}\;\;&h(B)=\frac{1+g_{0}}{G_{0}-2g_{0}}e^{\frac{G_{0}+2}{2k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}(B-B_{0})(B+B_{0}-2k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}})}.\end{split}$
(11)

The form of the equation for $B$ in Eq.(9) means that the magnitude of $B$ can
only ever decrease from the initial value of $B_{0}$. However, using Eq.(11),
we can determine the values of the initial conditions, and
$\\{k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}},k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\\}$
which allow $B$ to reach a value of zero. Due to the numerator of $g(B)$ in
Eq.(11), the Green sub-populations have the potential to be completely
depleted if the minimum value of $h(B)\leq 1/G_{0}$. As the argument of $h(B)$
in Eq.(11) is a positive quadratic, with minimum value
$B=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$, the minimum value of $h(B)$ is
$h(k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}})$. Thus, the condition between
the initial conditions and
$\\{k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}},k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\\}$
which enables $B\rightarrow 0$ is given by

$\displaystyle\begin{split}&h(k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}})>\frac{1}{G_{0}},\\\
\Rightarrow\;\;&k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}>B_{0}+\rho(k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}})-\sqrt{\rho(k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}})[\rho(k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}})+2B_{0}]},\end{split}$
(12)

where

$\rho(k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}})=\frac{k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}}{G_{0}+2}\ln\frac{G_{0}(1+g_{0})}{G_{0}-2g_{0}}.$
(13)

Thus, given initial conditions $G_{0}$, $g_{0}<G_{0}/2$ and
$B_{0}>k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$, the values of
$k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}$ and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$ which satisfy the requirement
given by Eq.(12) will ensure that $B\rightarrow 0$ after sufficient time.

Figure 5: Comparison of Eq.(12) with corresponding numerical outputs of final
$B$-values of Eqs.(4,6) with
$k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{L}}}_{\textrm{\tiny{B}}}$,
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$.
Initial conditions applied are $B_{0}=R_{0}=1.5$, $g_{0}=\gamma_{0}=1$,
$\Gamma_{0}=3$, resulting in system becoming Eq.(9) — i.e. $B=R$ and
$g=\gamma$. The black contour gives the case of Eq.(12) with an equals sign
which is the boundary between zero and non-zero $B$ values at the final time.

Figure 5 reveals very good agreement of Eq.(12) with corresponding numerical
outputs of final $B$-values stemming from Eqs.(4,6) with
$k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{L}}}_{\textrm{\tiny{B}}}$,
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$,
and matching initial conditions for $B,R$ and $g,\gamma$. Eq.(12) allows for
effective interrogation of the impact of changes to parameters and initial
conditions on the form of the contour which delineates between Green
contributing, or choosing to withdraw. For instance, if we mandate that
$|2B_{0}/\rho(k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}})|<1$, then performing
a binomial expansion on the square root of Eq.(12) obtains

$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}>\frac{B^{2}_{0}(G_{0}+2)}{2k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}}\ln\frac{G_{0}(1+g_{0})}{G_{0}-2g_{0}}}+O\left(\frac{B^{3}_{0}}{(k^{\textrm{\tiny{L}}}_{\textrm{\tiny{R}}})^{2}}\right).$
(14)

Eq.(14) reveals the relationships between carrying capacities and lethalities
for peer forces pursuing mutual annihilation with support from from the non-
combatant population — a situation to be avoided.

###  3.3 Contributing case: parameter sweeps

We offer phase plots of final population values for the Green-contributing
scenario in Figure 6, over the range
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\in(0,3.5)$ (Blue’s transfer
coefficient) and $k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\in(0.9,1.4)$
(Red’s capacity to maintain support from $g$). A narrower parameter range than
those considered in Figure 4 is given due to the contours being significantly
more complicated in the contributor scenario. The top-left panel for $R-B$
values at $t=50$ generally has three distinct regions: top-left diagonal of
Red victories; bottom-right diagonal of Blue victories; and middle-diagonal
consisting of a complex region of mixed outcomes separated by a highly
deformed spiral.

Figure 6: Parameter scans for Model II. The solutions are produced for
$B(t=0)=R(t=0)=20$, $g(t=0)=\gamma(t=0)=\Gamma(t=0)=10$, and
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1.2$.
Colours indicate the final population values (at $t=50$) over a range of
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}$ and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$, with blue and red signifying $B$
and $R$ victories in the top left panel. The white regions in the top-left
panel separate the Red and Blue winning regions.

In the top-left and bottom-right diagonal regions, trajectories generally
follow the pattern of Green contributors initially retreating to the neutral
$\Gamma$ population due to Red and Blue being over their respective carrying
capacity values. Thus Red and Blue initially damage each other at the same
rate as their starting populations are equal. At some point in time in either
of these regions, either the Red or Blue population becomes lower than their
carrying capacity value. Once this happens, the following outcomes are
possible:

  * •

In the upper-left diagonal region with
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}>k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$,
above the white contours, Red’s capacity to maintain support from Green is
sufficiently greater than Blue’s, and Blue’s transfer coefficient
($k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}$) is sufficiently less than Red’s.
This leads to a sharp rise in Green population $g$ who contribute to Red’s
effort, which in turn leads to a rise in Red’s population which is enough to
secure victory against Blue before any support for Blue can take effect.

  * •

For
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}<k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$
in the same upper-left diagonal region, Blue’s population decreases more
rapidly, leading to $\gamma$ increasing in population relatively early in the
engagement. However, because Blue’s transfer coefficient
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}$ is so small, the majority of the
corresponding contributors from Green never get the chance to actually
contribute. Hence, Red emerge victorious as they are much more effective at
transferring their Green contributors into the engagement after they start
receiving Green support.

  * •

In the lower-right diagonal region with
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}>k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$,
below the white contours, Red’s capacity to obtain contribution from Green is
greater than Blue’s. Thus Red’s contributor population from Green begins
increasing before Blue’s. Nevertheless in this region Blue’s transfer
coefficient ($k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}})$ is sufficiently
large enough that when Blue finally does obtain contributions from Green,
Blue’s ability to transfer said contributors into the engagement is enough to
counter Red’s efforts and win the engagement — even though Blue’s action was
slower compared to Red’s.

  * •

For
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}<k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$
in the lower-right diagonal region, as with the equivalent upper-left diagonal
region, Blue reaches the capacity to obtain contributions from Green before
Red, meaning that $\gamma$ increases in population before $g$. Blue’s transfer
coefficient $k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}$ is also sufficiently
large that the Green contributors are transferred effectively to swiftly
ensure Blue’s victory.

Generally in these regions above and below the white contours, the initial and
final total Green populations do not differ substantially from each other — as
seen by the top-right panel in Figure 6 for final $\Gamma$ population which
does not differ greatly from $G_{0}$ in these regions.

Figure 7: Number of periods in $R$ — signified by different colours ranging
from $0$ (black) to $7$ (yellow) periods. Numerical outputs are produced for
$B(t=0)=R(t=0)=20$, $g(t=0)=\gamma(t=0)=\Gamma(t=0)=10$. Left panel:
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1.2$.
Middle panel:
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1.2$.
Right panel:
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1.2$.

The trajectories of the final region inside the white contours, centred at the
point
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}=k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1.2$,
the similarity of the parameter values means the trajectories for the
engagements between Red and Blue display a complex back-and-forth dynamic as
the population of each force increases and decreases a number of times. Figure
7 displays the number of periods experienced by Red’s trajectory as a function
of input parameters, with the middle panel corresponding to the parameter
inputs used in Figure 6. Focusing on the middle panel in Figure 7, the region
bounded by the white contours exhibits significant oscillatory behaviour —
centred at
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}=k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}=1.2$.
Mirroring the top two panels of Figure 3, the oscillations of the Red and Blue
trajectories are centred around each force’s carrying capacity value. There is
a stark contrast when comparing the final total Green population(s) in this
region with the regions above and below the white contours. The oscillations
generally serve to prolong the engagement, as seen by the difference between
the time-scales of the top and bottom panels of Figure 3. Thus the final total
Green populations in the oscillatory region is notably much less than the
surrounding regions as the prolonged time period means that a greater
proportion of the Green population contributes to the conflict.

###  3.4 Contributing case: analytic considerations

By applying approximations in this instance of the model it is possible to
gain greater understanding of the oscillatory behaviour displayed by the Red
and Blue trajectories. As with the previous case in section 3.2 where the
near-peer scenario was examined in depth, to allow analytical amenability we
assign
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}=k^{\textrm{\tiny{T}}}_{\textrm{\tiny{R}}}$
and
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{B}}}$
in Eq.(8), and consider the case $\gamma=g$ and $B=R$. Doing so to we obtain
the reduced system:

$\displaystyle\dot{B}=-B+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}g,\;\;\dot{g}=g\Gamma
B\left(1-\frac{B}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}\right)-k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}g,\;\;\dot{\Gamma}=-2g\Gamma
B\left(1-\frac{B}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}\right).$ (15)

Further assuming that $B$ exhibits small oscillations around the value
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$, behaviour that was exhibited in
the top two panels of Figure 3. Hence,

$B(t)=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}+\epsilon(t),$ (16)

where we assume that the time dependent variable in Eq. (16) is small, i.e.
$\epsilon^{2}\approx 0$. With this assumption, the system in Eq. (15) becomes,

$\dot{\epsilon}=-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}-\epsilon+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}g,\;\;\dot{g}=-g\Gamma\epsilon-k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}g,\;\;\dot{\Gamma}=2g\Gamma\epsilon,$
(17)

where any quadratic terms in $\epsilon$ are approximately zero. The
multiplicative $\epsilon$ terms on the right hand side of Eq. (17) highlight
their role in determining the direction of flow between $g$ and $\Gamma$. That
is, if $\epsilon<0$, Blue can maintain support from Green, and population will
flow from $\Gamma$ to $g$. Contrastingly, if $\epsilon>0$, Blue cannot
maintain support and contributions from Green, and population will flow from
$g$ to $\Gamma$.

Summing all the terms in Eq. (17), and assuming that the initial condition for
$\epsilon$ is zero (i.e. $\epsilon_{0}=0$), we obtain

$\displaystyle\begin{split}&\dot{\epsilon}+\dot{g}+\frac{1}{2}\dot{\Gamma}=-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}-\epsilon,\\\
\Rightarrow\;\;&\epsilon+g+\frac{1}{2}\Gamma=g_{0}+\frac{1}{2}\Gamma_{0}-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}t-\int^{t}_{0}d\tau\epsilon(\tau).\end{split}$
(18)

In the dynamical region of interest of oscillatory $\epsilon$, we can assume
the term $\int^{t}_{0}d\tau\epsilon(\tau)$ contributes very little to Eq.(18)
as the trajectory of $\epsilon$ will not have strayed significantly negatively
away from zero. Thus at $t=t_{f}$, where

$t_{f}=\frac{1}{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}\left(g_{0}+\frac{\Gamma_{0}}{2}\right),$
(19)

the right hand side of the final expression of Eq. (18) approximately equals
zero, giving a natural time-scale, past which the oscillatory dynamics will
cease, as either both $g$ and $\Gamma$ are zero, or $\epsilon$ is
significantly below zero.

Dividing $\dot{g}$ by $\dot{\Gamma}$ in Eq. (17) reveals

$\displaystyle\begin{split}2\frac{dg}{d\Gamma}=-1-\frac{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}}{\Gamma\epsilon},\\\
\Rightarrow-2\epsilon(d\epsilon+k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}dt+\epsilon
dt)=-k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\frac{\ d\Gamma}{\Gamma},\\\
\Rightarrow\epsilon^{2}+2\int^{t}_{0}d\tau\left[k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\epsilon(\tau)+\epsilon^{2}(\tau)\right]=k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\ln\frac{\Gamma}{\Gamma_{0}},\\\
\Rightarrow\Gamma=\Gamma_{0}\exp\left\\{\frac{\epsilon^{2}+2\int^{t}_{0}d\tau\left[k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\epsilon(\tau)+\epsilon^{2}(\tau)\right]}{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}}\right\\}\end{split}$
(20)

where Eq. (18) has been applied to proceed from the first line of Eq. (20) to
the second. Additionally, the expression $\int^{t}_{0}d\tau\epsilon(\tau)$,
and the terms nonlinear in $\epsilon$ in Eq.(20) have been kept, with a plan
to deal with them appropriately at the conclusion of this section. Eq. (20)
gives the trajectory of $\Gamma$ entirely in terms of $\epsilon$. In order to
utilise this expression, we re-cast the equation for $g$ in Eq. (17) as

$e^{-k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}t}\frac{d}{dt}ge^{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}t}=\dot{g}+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}g=-\frac{1}{2}\dot{\Gamma}.$
(21)

Performing a similar operation to the equation for $\epsilon$ in Eq. (17) —
i.e. multiplying the entire expression by
$e^{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}t}$, and then differentiating
with respect to time — we obtain

$\displaystyle\begin{split}\ddot{\epsilon}+(k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}+1)\dot{\epsilon}+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\epsilon=-\frac{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}}{2}\dot{\Gamma}-k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\\\
\Rightarrow\dot{\epsilon}+(k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}+1)\epsilon+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\int^{t}_{0}d\tau\epsilon(\tau)+\frac{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}}{2}\Gamma=k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\left(g_{0}+\frac{\Gamma_{0}}{2}-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}t\right)-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}},\end{split}$
(22)

with initial condition
$\dot{\epsilon}_{0}=-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}g_{0}$.
As the expression for $\Gamma$ in terms of $\epsilon$ in Eq. (20) contains an
integral of $\epsilon$ with respect to time, the expression in Eq. (22) gives
an integro-differential equation for the trajectory of $\epsilon$. To change
all but the argument of the exponential of Eq. (22) into an ODE, we redefine
the time dependent variables as follows,

$\alpha=\int^{t}_{0}d\tau\epsilon(\tau),\;\;\dot{\alpha}=\epsilon,\;\;\ddot{\alpha}=\dot{\epsilon},$
(23)

with initial conditions $\alpha_{0}=\dot{\alpha}_{0}=0$,
$\ddot{\alpha}_{0}=-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}g_{0}.$
Doing so, Eq. (22) becomes:

$\ddot{\alpha}+\frac{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\Gamma_{0}}{2}\exp\left\\{\frac{\dot{\alpha}^{2}+2\left[k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\alpha+\int^{t}_{0}d\tau\dot{\alpha}^{2}(\tau)\right]}{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}}\right\\}+(k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}+1)\dot{\alpha}+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\alpha=k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\left(g_{0}+\frac{\Gamma_{0}}{2}-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}t\right)-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}.$
(24)

The expression for the trajectory of $\epsilon$ — i.e. $\dot{\alpha}$ — in Eq.
(24) produces exactly the trajectories for $\epsilon$ numerically generated
from the system in Eq. (17). Nevertheless, being nonlinear itself Eq. (24)
still requires numerical integration.

In order to gain analytical insight, we first approximate Eq. (24) into an ODE
by ignoring the integral term in the exponent to obtain,

$\ddot{\alpha}+(k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}+1)\dot{\alpha}+k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\alpha+\frac{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\Gamma_{0}}{2}\exp\left(\frac{\dot{\alpha}^{2}+2k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\alpha}{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}}\right)=k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}\left(g_{0}+\frac{\Gamma_{0}}{2}-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}t\right)-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}.$
(25)

which is of the form of a forced autonomous equation — refer to section 2.2 of
(Polyanin and Zaitsev 2003) for more details. Next we assume that the
quadratic term of $\dot{\alpha}$ in exponential of Eq.(25) equals zero (i.e.
$\dot{\alpha}^{2}\equiv\epsilon^{2}\approx 0$). Additionally, we linearise the
remaining term in the exponential via,

$\exp\left(\frac{2k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}}\alpha\right)\approx
1+\frac{2k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}{k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}}\alpha.$
(26)

Doing so, Eq. (25) becomes,

$\displaystyle\begin{split}\ddot{\alpha}+(k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}+1)\dot{\alpha}+(k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}+k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\Gamma_{0})\alpha=k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}(g_{0}-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}t)-k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}.\end{split}$
(27)

In the left and right panels of Figure 8 we plot time-dependent trajectories
for $B$ and $\epsilon$, and the $\Gamma$ vs $B$ and $\epsilon$ phase-space
trajectories for the system Eq. (15), and its various approximations given by
Eqs. (17), (25) and (27). Importantly, the solutions for $\epsilon$ are
shifted by $k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$ for comparison as per
Eq. (16). On the left panel we see that all four time-dependent trajectories
initially produce very similar profiles. Nevertheless (and perhaps
unsurprisingly) the trajectory corresponding to linearised Eq.(27) deviates
significantly from the other three profiles, unable to reproduce the
nonlinearity clearly being demonstrated for the entirety of the time-scale
Eq.(19).

Figure 8: Solutions for $B$ and $\epsilon$ (left panel), and phase portraits
(right panel) for the systems Eqs. 15, 17, and Eqs. 25 and 27, are shown by
black, red, blue and magenta curves, respectively. The solution is produced
for $B_{0}=B(t=0)=1$, $g_{0}=g(t=0)=3$, $\Gamma_{0}=\Gamma(t=0)=20$, and
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}=k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}=1$.
The solutions for $\epsilon$ are shifted by
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$ for comparison (see Eq. 16).

Thus, although the linearity of Eq.(27) clearly limits its predictive ability
over the entire time-scale, it is nonetheless immediately utilisable as the
general solution is a sum of a linear combinations of exponents and a linear
function. Importantly, oscillatory solutions exist if the discriminant of the
characteristic equation is negative, i.e.

$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}<1+2\sqrt{k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\Gamma_{0}},$
(28)

and the oscillation period is then given as

$T=\frac{2\pi}{\omega}=\frac{4\pi}{\sqrt{|(k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}-1)^{2}-4k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\Gamma_{0}|}}.$
(29)

Therefore, the number of oscillation periods in the solution can be
approximated as the ratio of the timescale, given in Eq. (19), to the
oscillation period as:

$N\simeq\frac{t_{f}}{T}=\frac{\left(2g_{0}+\Gamma_{0}\right)\sqrt{|(k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}-1)^{2}-4k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}\Gamma_{0}|}}{8\pi
k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}}.$ (30)

Hence the combination of Eq.(30) and (29) enables an approximate understanding
of the effect that changes in input parameters and initial conditions have on
the applicable regions, and number of oscillations occurring within the the
time-scale.

Figure 9: Number of periods calculated from Eq. 30 (white numbered contours)
in comparison to the number of periods in the solution of the system Eq. 15 at
final time, given by the coloured regions ranging from black ($0$ periods) to
yellow ($5$ periods). The initial conditions are as in Fig. 8.

In Figure 9 we compare the oscillatory regions, and number of periods,
generated from both the original system in Eq.(15) — represented by the
coloured contours — and the fully linearised approximations given by equations
(28) and (30) — represented by the white numbered contours. The contours
generated by Eq.(30) do display similar macroscopic properties exhibited by
the system in Eq.(15). Notably, that an increase in oscillations is
approximately inversely proportional to the carrying capacity value
$k^{\textrm{\tiny{C}}}_{\textrm{\tiny{R}}}$. Additionally, we also see that as
the transfer rate $k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}$ decreases, this
generally leads to an increase in oscillation count — though from the contours
generated by Eq.(15), we see that a further decrease in
$k^{\textrm{\tiny{T}}}_{\textrm{\tiny{B}}}$ leads to the oscillation count
decreasing, something not reflected in the contours of the linear system.

##  4 Discussion and future work

In this work we have extended the classical adversarial Blue vs Red model of
combat attrition to include a third neutral population who are given the
opportunity to become supporters and add to the combat effectiveness of Red
and Blue, or contributors who seek to actively join Red and Blue in their
engagement. Blue or Red victory in the supporter model was most sensitive to
changes of the carrying capacity. Notably, if the carrying capacity value was
less than the corresponding force’s initial population, the likely outcome was
defeat for that force in the absence of a change of strategy. We can see such
sensitivity to carrying capacity in the aftermath of the 2003 Iraq war. The
failure to prevent looting in the immediate days, followed by the later de-
Ba’athification triggering the insurgency (Hosmer 2007) may be interpreted as
a collapse in carrying capacity for US forces within its potential supporting
population. The subsequent 2007 Surge under the ‘Petraeus doctrine’ (Petraeus
2009) represented both a numerical increase in forces to defeat insurgents and
the provision of security to local populations to enhance internal governance,
effectively to increase the carrying capacity in potential supportive
communities.

The contributor model revealed regions of clear Blue and Red victories which
surrounded a highly volatile region where outcome prediction were largely
impossible. This volatility was due to both Blue and Red, and their
corresponding Green supporters, being relatively evenly matched, leading to a
long protracted conflict. During this conflict Blue and Red population values
oscillated about their carrying capacity values, fueled by the continual
contribution of Green population into both Blue and Red combatants, until
finally being depleted. This scenario is reminiscent of the early stages of
the 2014 conflict in Eastern Ukraine, where near-peer Ukrainian forces and
pro-Russian separatists received an influx of contributors from surrounding
regions (Kofman et al. 2017).

Under sufficiently symmetrised assumptions, corresponding to a worst-case
stalemate scenario, for the supporter variant we related the dynamic
variables, initial conditions and the input parameters. This relationship
offered a computationally inexpensive means to determine if the engagement
between Red and Blue would end due to Green rescinding support completely in a
near-peer setting. For the contributor model, under similar assumptions which
sufficiently symmetrised defining equations, we were able to analyse the
protracted oscillatory regime where a near-peer Red and Blue fluctuated around
the carrying capacity value. Notably, we were able to derive a closed form
expression of the expected time-scale of the oscillations, and a non-linear
ODE describing the oscillations, both of which performed quite well when
compared numerically with the original system.

As the original intent of this work was to devise a relatively simple model
which offered a means to explore the mutual influence between combatants and
non-combatants in warfare, we shall devote some effort in interpreting these
model behaviours through this lens; firstly by discussing possible meanings
and implications for each force’s carrying capacity value. Carrying capacity
was the means we introduced each Green sub-population’s growth and decline of
support. If either Blue or Red grew too big, their potential Green followers
retract support. This is to be interpreted as the Green population perceiving
either force as acting as a military occupation, and the negative consequences
which stem from that. Consequences which may be as devastating as collateral
damage, or as the negative perception held by a local population that either
Red or Blue are exerting too much influence on local decision-making. Hence, a
higher carrying capacity value for either Red or Blue equates to that force
being more disciplined so as not to cause collateral damage, and willing to
work with a local population, allowing for self-determination and not denying
their ability for meaningful decision-making. Such a force would possess
superior training, both military and cultural, and be steeped in doctrine
which does not shy away from distributed and/or collaborative decision-making.

Bringing this back to behaviours witnessed by the model, we see that if it is
Red and/or Blue’s intent to ensure a decisive victory, then figures 4 and 6
clearly show that the method to ensure this through possessing a superior
carrying capacity (training, doctrine, decision-making, trust, etc.) for
either supporter or contributor model variant, as this parameter is the most
sensitive to change in model outcomes once approximate parity is reached. In
the supporter model, with ensured conservation of Green, it would be sensible
for either Blue or Red to aim for advantageous regions if there is approximate
parity in carrying capacities — equivalent to the band of Red victory
surrounded by Blue victories on either side in figure 4. If however there is
possibility for Green involvement in the engagement, similar to the
contributor model variant, then it is paramount to avoid approximate parity in
carrying capacity values as this leads to protracted engagements with
potentially disastrous consequences on the Green population as witnessed by
the depleted population values in figure 6.

For future applications, an immediate generalisation would include a mixed
variant of our model, with one force being supported and one being contributed
to, similar to Deitchman’s (1962) asymmetric guerrilla warfare model for the
original Lanchester equations. It is our intention nevertheless to couple this
model with a system of networked decision-makers via the Kuramoto-Sakaguchi
model considered in (Zuparic et al. 2021; Kalloniatis et al. 2020). Ideally
such a coupled model would possess the complexities exhibited in the model
considered in this work, with the advantage of having the ability to more
effectively unpack the interpretation of the carrying capacity term into its
implied combat and decision-making meanings.

## Acknowledgements

This research was a collaboration between the Commonwealth of Australia
(represented by the Defence Science and Technology Group) and Deakin
University through a Defence Science Partnerships agreement, under the
auspices of the Modelling Complex Warfighting initiative.

## References

  * Atkinson et al. (2012) Atkinson M., Gutfraind A. and Kress M. (2012) When do armed revolts succeed: lessons from Lanchester theory. Journal of the Operational Research Society. 63(10):1363–73.
  * Arató (2003) Arató M. (2003) A famous nonlinear stochastic equation (Lotka-Volterra model with diffusion). Mathematical and Computer Modelling 38(7–9):709–26.
  * Deitchman (1962) Deitchman S. (1962) A Lanchester Model of Guerrilla Warfare. Operations Research. 10(6):818–27.
  * Dorn and Varey (2009) Dorn W. and Varey M. (2009) Fatally flawed: The rise and demise of the “three-block war” concept in Canada. Canadian Military Journal. 10(1):38–45.
  * Drapeau (2008) Drapeau, M., Hurley, P. and Armstrong, R. (2008) So many zebras, so little time: Ecological Models and Counterinsurgency Operations. Defense Horizons 62: 1–8.
  * González and Villena (2011) González E. and Villena M. (2011) Spatial Lanchester models European Journal of Operational Research. 210(3):706–15.
  * González-Avella et al. (2014) González-Avella J., Cosenza M. and San Miguel M. (2014) Localized coherence in two interacting populations of social agents. Physica A. 399:24–30.
  * Hening and Nguyen (2018) Hening A. and Nguyen D. (2018) Persistence in stochastic Lotka-Volterra food chains with intraspecific competition. Bulletin of Mathematical Biology 80(10):2527–60.
  * Hening et al. (2018) Hening A., Nguyen D. and Yin G. (2018) Stochastic population growth in spatially heterogeneous environments: the density-dependent case. Journal of Mathematical Biology 76(3):697–754.
  * Holling (1966) Holling C. (1966) The functional response of invertebrate predators to prey density. Memoirs of the Entomological Society of Canada. 98(S48):5–86.
  * Hosmer (2007) Hosmer S. (2007) Why the Iraqi resistance to the coalition invasion was so weak. Santa Monica CA: RAND Corporation, 2007.
  * Hung et al. (2005) Hung, C., Yang G., Deng P., Tang T., Lan S. and Chu P. (2005) Fitting Lanchester’s square law to the Ardennes Campaign. Journal of the Operational Research Society 56(8):942–6.
  * Johnson (2008) Johnson D. and Madin J. (2008) Population Models and Counterinsurgency Strategies. In Natural Security: A Darwinian Approach to a Dangerous World, Eds. R. D. Sagarin T. Taylor, University of California Press, 159.
  * Kalloniatis et al. (2021) Kalloniatis A., Hoek K., Zuparic M. and Brede M. (2021) Optimising structure in a networked Lanchester model for Fires and Manoeuvre in Warfare. Journal of the Operational Research Society. 72(8):1863–78.
  * Kalloniatis et al. (2020) Kalloniatis A., McLennan-Smith T. and Roberts D. (2020) Modelling distributed decision-making in Command and Control using stochastic network synchronisation. European Journal of Operational Research. 284(2):588–603.
  * Keane (2011) Keane T. (2011) Combat modelling with partial differential equations. Applied Mathematical Modelling. 35(6):2723–35.
  * Kilcullen (2020) Kilcullen D. (2020) The Dragons and the Snakes: How the rest learned to fight the west. Melbourne: Scribe.
  * Kofman et al. (2017) Kofman M., Migacheva K., Nichiporuk B., Radin A., Tkacheva O. and Oberholtzer J. (2017) Lessons from Russia’s operations in Crimea and Eastern Ukraine. Santa Monica CA: RAND Corporation, 2017.
  * Kress and Szechtman (2009) Kress M. and Szechtman R. (2009) Why defeating insurgencies is hard: The effect of intelligence in counterinsurgency operations——A best-case scenario. Operations Research. 57(3):578–85.
  * Kress et al. (2018) Kress M., Lin K. and MacKay N. (2018) The attrition dynamics of multilateral war. Operations Research. 66(4):950–6.
  * Kress and MacKay (2014) Kress M. and MacKay N. (2014) Bits or shots in combat? The generalized Deitchman model of guerrilla warfare. Operations Research Letters. 42(1):102–8.
  * Lanchester (1916) Lanchester F. (1916) Aircraft in warfare: Dawn of the fourth arm. London: Constable Limited.
  * Lin et al. (2014) Lin K. and MacKay N. (2014) The optimal policy for the one-against-many heterogeneous Lanchester model. Operations Research Letters. 42(6-7):473–7.
  * Lotka (1925) Lotka A. (1925) Elements of physical biology. Maryland: Williams and Wilkins.
  * Lucas and Turkes (2004) Lucas T. and Turkes T. (2004) Fitting Lanchester equations to the battles of Kursk and Ardennes. Naval Research Logistics. 51(1):95–116.
  * MacKay (2006) MacKay N. (2006) Lanchester combat models. Mathematics Today. 42(5):170–3.
  * MacKay (2009) MacKay N. (2009) Lanchester models for mixed forces with semi-dynamical target allocation. Journal of the Operational Research Society. 60(10):1421–7.
  * MacKay (2012) MacKay N. (2012) Response to ’Some comments on MacKay’s Lanchester mixed forces model’. Journal of the Operational Research Society. 63(11):1633.
  * MacKay (2015) MacKay N. (2015) When Lanchester Met Richardson, the Outcome Was Stalemate: A Parable for Mathematical Models of Insurgency. In: Forder R, editor. OR, Defence and Security. London: Palgrave Macmillan. p. 124–47.
  * McLennan-Smith et al. (2019) McLennan-Smith T., Kalloniatis A., Sidhu H., Jovanoski Z., Watt S. and Towers I. (2019) Exploiting ecological non-trophic models in representations of warfare. Elsawah S, editor. MODSIM2019, 23rd International Congress on Modelling and Simulation. Modelling and Simulation Society of Australia and New Zealand, December. pp. 242–-8.
  * McLennan-Smith et al. (2020) McLennan-Smith T., Kalloniatis A., Jovanoski Z., Sidhu H., Roberts, D., Watt S. and Towers I. A mathematical model of humanitarian aid agencies in attritional conflict environments. Operations Research. 69(6): 1696–-1714.
  * Metz and Millen (2004) Metz S. and Millen R. (2004) Insurgency and counter-insurgency in the 21st Century: Reconceptualizing threat and response. Strategic Studies Institute, US Army War College
  * Osipov (1995) Osipov M. (1995) The Influence of the Numerical Strength of Engaged Forces on Their Casualties (translated by Helmbold R and Rehm A). In: Bracken J, Kress M and Rosenthal R, editors. Warfare Modelling. Massachussets: John Wiley and Sons. p. 289–344.
  * Petraeus (2009) Petraeus D. (2009) The Petraeus Doctrine: The Field Manual on Counterinsurgency Operations Joint Chiefs of Staff Joint Publication 3-24.
  * Polyanin and Zaitsev (2003) Polyanin A. and Zaitsev V. (2003) Handbook of Exact Solutions for Ordinary Differential Equations (2nd Edition). Boca Raton: Chapman and Hall/CRC.
  * Protopopescu et al. (1989) Protopopescu V., Santoro R. and Dockery J. (1989) Combat modeling with partial differential equations. European Journal of Operational Research. 38(2):178–83.
  * Richards (1959) Richards F. (1959) A flexible flexible growth function for empirical use. Journal of Experimental Botany 10(2):290–301.
  * Royama (1971) Royama T. (1971) A comparative study of models for predation and parasitism. Population Ecology 13:1–91.
  * Schaffer (1968) Schaffer M. (1968) Lanchester Models of Guerrilla Engagements. Operations Research. 16(3):457–88.
  * Syms and Solymar (2015) Syms R. and Solymar L. (2015) A dynamic competition model of regime change. Journal of the Operational Research Society. 66(11):1939–47.
  * Taylor (1974) Taylor J. (1974) Solving Lanchester-type equations for “modern warfare” with variable coefficients. Operations Research. 22(4):756–70.
  * Taylor and Parry (1975) Taylor J. and Samuel P. (1975) Force-ratio considerations for some Lanchester-type models of warfare. Operations Research 23(3):522–33.
  * Taylor and Brown (1976) Taylor J. and Brown G. (1976) Canonical methods in the solution of variable-coefficient Lanchester-type equations of modern warfare. Operations Research. 24(1):44–69.
  * Taylor (1977) Taylor J. (1977) Determining the class of payoffs that yield force-level-independent optimal fire-support strategies. Operations Research. 25(3):506–16.
  * Taylor (1979) Taylor J. (1979) Optimal commitment of forces in some Lanchester-type combat models. Operations Research. 27(1):96–114.
  * Verhulst (1838) Verhulst P. (1838) Notice sur la loi que la population suite dans son accroissement. Correspondence Mathematique 10:113–21.
  * Volterra (1928) Volterra V. (1928) Variations and fluctuations of the number of individuals in animal species living together (translated by Wells M). ICES Journal of Marine Science. 3(1):3–51.
  * Williamson (2011) Williamson, J. (2011) Using humanitarian aid to win hearts and minds’: a costly failure? International Review of the Red Cross 93(884):1035.
  * Zimmerman (2016) Zimmerman D.(MAJ) (2016) Understanding the Threat Ecosystem: A Concept for Intelligence Support to Special Warfare. School of Advanced Military Studies United States Army Command and General Staff College Fort Leavenworth, Kansas.
  * Zuparic et al. (2021) Zuparic M., Angelova M., Zhu Y. and Kalloniatis A. (2021) Adversarial decision strategies in multiple network phased oscillators: the Blue-Green-Red Kuramoto-Sakaguchi model. Communications in Nonlinear Science and Numerical Simulation. 95:105642.