A size-structured matrix model to simulate dynamics of marine community size spectrum.

Several types of size-based models have been developed to model the size spectra of marine communities, in which abundance scales strongly with body size, using continuous differential equations. In this study, we develop a size-structured matrix model, which can be seen as a discretization of the Mckendrick-von Foerster equation, to simulate the dynamics of marine communities. The developed model uses a series of simple body size power functions to describe the basic processes of predator-prey interactions, reproduction, metabolism, and non-predation mortality based on the principle of mass balance. Linear size spectra with slopes of approximately -1 are successfully reproduced by this model. Several examples of numerical simulations are provided to demonstrate the model's behavior. Size spectra with cut-offs and/or waves are also found for certain parameter values. The matrix model is flexible and can be freely manipulated by users to answer different questions and is executed over a shorter numerical calculation running time than continuous models.

Introduction

Since the pioneering work by Sheldon et al. [1], the marine community size spectrum, i.e., the distribution of abundance or biomass against individuals’ body size on a log–log scale, has attracted significant interest from the scientific community. Classifying individuals over the size range from bacteria to whales into different size groups regardless of taxonomy results in a linear size spectrum with a slope of approximately –1 in unexploited communities [1, 2]. The slopes and intercepts of these size spectra have been proposed as useful indicators in revealing the nature of ecosystems and the effects of anthropogenic perturbation on the structure of marine communities [3].

Several types of size-based models have been developed to reproduce linear size spectra, explain their remarkable regularities, and explore the related ecosystem dynamics [4-8]. Such models are also widely used to investigate the impacts of environmental changes [9, 10] and fishing [11-13] on marine communities. Such models have been used to suggest that different types of behaviors can theoretically emerge at various size spectra. Using a jump-growth model in which predation is represented as a Markov process, Datta et al. [14] demonstrated that the size spectrum can switch from a linear to a wave-like shape when predators feed over a relatively narrow size range of prey.

In these studies, the balance between mortality, growth, reproduction, metabolism, etc. at all body sizes is formalized by using the Mckendrick-von Foerster partial differential equation [15]. Partial differential equation models (PDE models) enable the incorporation of detailed processes with various parametric functions. Although the R package “mizer” developed by Scott et al. [8] has been widely used for implementing size-based models, it may not always be easy or convenient for users to couple processes such as environmental stochasticity to PDE models in obtaining fits to the real situations of specific marine communities [16]. Matrix models can be more easily manipulated than PDE models in modeling such processes because they use simple discrete states [16]. Matrix algorithms can also be executed within shorter running times than continuous models with numerous loop iterations. Because real data are always recorded in discrete form, matrix models can be calibrated and validated with real data points.

In marine ecosystems, body size correlates with various biological and physiological processes; for example, the vital rates at the individual level, including encounter rate, swimming speed, reproduction rate, metabolic rate, and natural mortality rate, all ubiquitously scale allometrically with body size [7, 17–20]. Such processes are considered to be the main drivers that affect the dynamics of abundance and, because they can be described by simple power law functions of body size, they can be easily incorporated into a matrix model. Models such as the Leslie matrix model [21] have been used in population ecology to predict population growth and the age composition of the population. These models can provide clues for the construction of a size-structured matrix model for marine communities. Motivated by this background, the main goal of this study was to develop a simple and widely applicable size-structured matrix model that can easily be modified and extended to suit different applications in flexibly investigating community- and population-level dynamics.

In the following, we first present matrices corresponding to the basic processes of community dynamics to construct a unified matrix model. We then provide several numerical simulation examples in which the values of the model parameters are varied separately in a range around their reference values. The effects of fishing on stationary size spectra are then predicted under various fishing scenarios. Finally, we discuss the model structure and behavior of the observed size spectra and provide some modifications of the model for future applications.

Materials and methods

Size-structured community dynamics

Consider a size-structured community (e.g., marine community) in which larger organisms prey on smaller organisms (in this paper, we do not discriminate the identities of species). We begin construction of the matrix model by discretizing individuals with different body masses into k equally-spaced logarithmic size classes. The community structure is represented as a column vector n(t) = (n1(t),n2(t),…,n(t))T, where the superscript T denotes vector or matrix transposition. Each value in the vector represents the number of individuals per unit volume measured in ind./m-3 within the class at time t (for brevity, in the following we omit writing (t) for the individual size classes). Let m = (m1,m2,…,m)T be the column vector describing the body masses of individuals in each size class. The body mass axis is discretized using a mass ratio of successive size classes δ: m = mδ, where i = 1,2,…,k. In our model, the community is divided into two groups: the background resources’ group, which has body masses m1 and produces the only input to the community; and the consumers’ group, which has body masses belonging to [m2,m].

The model, presented graphically in Fig 1, shows four processes that affect the number of individuals in each size class per unit time Δt. The dynamics of community size structure can be described in terms of these basic processes of predator–prey interaction, reproduction, metabolism, and non-predation mortality. As described below, the model presented here relies on empirical knowledge (e.g., allometric relationships) between body mass (m) and the rates or quantities of organism physiological phenomena (y) according to simple power functions y = am, where a is a coefficient and b is a scaling exponent. Our model is designed based on a mass balance in the community (i.e., no biomass is wasted). Because materials such as dead animals, fecal pellets, and metabolites from consumers constitute the assimilable inorganic and organic nutrients needed for phytoplankton growth, we assume that all egested materials (i.e., undigested materials, metabolites, carrion, etc.) are assimilated to form new resources through microbial degradation such that a simple modeling of the dynamics of resources is possible. To maintain the mass balance, an imaginary size class k+1 is introduced to remove the increased biomass in size class k induced by predation.

Fig 1

Graphic representation of the size-structured matrix model.

Nodes represent size classes; larger circles depict larger body mass. Background resources are assumed to occupy only the first size class and consumers to range from the second to the kth size class. Thin arrows indicate flow of matter and energy expressed in mass (E, C, F, W) or number of individuals (P, D, Y) between size classes through each process per unit time Δt, while thick arrows indicate transfer of individuals (G, B, L) between size classes owing to increase or decrease in body mass. Individuals can grow or move down no more than a single size class per unit time Δt. (a) Predator–prey interactions: in a predation process, several individuals P in prey size class j are eaten by individuals in predator size class i. Undigested materials E are egested and recycled to resources in size class 1. After feeding, a certain fraction of individuals G in predator size class i grow to the next size class i+1 owing to mass increase. For individuals in the largest size class k in this community, we assume that no predation occurs (i.e., P is zero for all i). Instead, we assume that all individuals growing from size class k into an imaginary larger and empty size class k+1 out of this community die and the carrion C is recycled to resources. (b) Reproduction: the first size class of the consumers’ group (i.e., size class 2) is represented as the starting point of consumers as eggs or larvae. Males and females reproduce at the same reproduction rate F, producing eggs into size class 2 (females) and excretive matters (semen) recycled to resources (males). As a consequence, a certain fraction of individuals B in parent size class i decrease their body mass and are transferred to the next-lower size class i–1. (c) Metabolism: individuals within size classes 2 to k produce metabolites that are recycled to resources at a rate W, with a certain fraction L transferred into the next-lower size class owing to loss in body mass. (d) Non-predation and fishing mortality: several individuals D die owing to natural (non-predation) mortality; meanwhile, their carrion is recycled to resources. Fishing removes some individuals Y from each size class i. These catches are consumed by humans and eventually reintegrated into the recycling process and back into the ecosystem as nutrients to resources. In our model, they are recycled in size class 1 while maintaining the mass balance.

Assuming that the basic processes in Fig 1 are independent of each other, the column vector n(t+Δt) representing the community structure at time t+Δt is given by where Δn, Δn, Δn, and Δn are column vectors (, , , and ) describing shifts in the number of individuals within each size class as a result of predator–prey interactions (P), reproduction (R), metabolism (Q), and non-predation, and fishing mortality (U), respectively.

The matrix model

Predator–prey interactions

In each predation event, we assume that predators have a preferred predator–prey mass ratio (PPMR) based on the rule “bigger ones eat smaller ones” [22]. A size selectivity matrix S by a predator in size class i on a prey in size class j composed of size preferences φ (i>j) and zeros is described as follows: where φ is given by a function that follows a lognormal distribution on the ratio of body masses of predator and prey, m and m, respectively [22]. Although resources are only placed at the first size class [m1] in our model, we assume that they also extend to an imaginary size range [m−,m0] (which is divided into logarithmically equally spaced s+1 size classes, and m0

Table 1

Parameter values used in the basic model simulation case.

Parameter	Definition	value	Source
γ	Factor of searching rate	600 m-3 yr-1 g-p	[23]
p	Exponent of searching rate	0.75	[13]
υ	Factor of preferred PPMR	100	[22]
Τ	Exponent of preferred PPMR	0	[22]
σ	Width of size preference	1.3	[22]
v	Factor of standard metabolic rate	2.4 yr-1 g1-q	[13]
q	Exponent of standard metabolic rate a	0.75	[24]
μ0	Factor of non-predation mortality	0.6 yr-1 g-d	[11]
d	Exponent of non-predation mortality	0.25	[19]
ks	Factor of carrying capacity a	0.0015	Modified from [4]
f	Exponent of carrying capacity b	1.15	Modified from [4]
c	Factor of reproduction rate c	10 yr-1 g1-r	Estimates from [6]
r	Exponent of reproduction rate d	0.75
α	Sex ratio e	0.5	Arbitrary
e	Assimilation efficiency f	0.6	[13]
q0	Catchability coefficient	0.1	Arbitrary
Δt	Time step g	0.001 yrs
δ	Mass ratio of successive size classes	1.5

aAdjusted from the factor of natural mortality rate in [4] to output a linear size spectrum.

bAdjusted from the exponent of natural mortality rate in [4] to output a linear size spectrum.

cEstimated from the fraction of the assimilated energy allocated to reproduction.

dWe use the value r = p = q = 0.75 assuming that the reproduction rate is proportional to the searching rate and the standard metabolic rate.

eThe proportion of females (sex ratio) is assumed to be constant across all size classes, i.e., α1 = α2 = ⋯ = α.

fIndividuals in each size class i are assumed to have the same assimilation efficiency, i.e., e1 = e2 = ⋯ = e.

gIt should be ensured that the probability of each individual growing up to its next size class can be no larger than one in each time step.

Suppose that the volume of water searched by a predator per unit time Δt is expressed as an allometric function of predator body mass γmΔt [25] and that the encounter rate of each predator individual to prey in size class j is also proportional to the number of individuals in prey size class j. Then, the total number of individuals in a prey size class j eaten by predators (∑P) per unit time Δt can be expressed as where N = diag(n), ST is the transpose of matrix S, and , describing the search rate of each predator individual in size class i. Note that the elements on the diagonal matrix N represent the number of prey individuals in each size class, whereas the entries of n represent the number of predator individuals.

On the other hand, the total biomass of individuals assimilated by predators in size class i with an assimilation efficiency e (e∑mP) per unit time Δt is calculated as where E = diag(e1,e2,…,e), whose entry e describes the assimilation efficiency of individuals in each predator size class i and M = diag(m) represents the prey body mass. Note that, in this case, the diagonal matrix N represents the number of predator individuals, while the vector n represents the number of prey individuals, which is different from the case of Eq (4). The first entry of the product ENVSMn becomes zero, which means that resources conceptually do not prey on others.

Increases in the body mass of predators in size class i owing to assimilation of prey can cause a certain number of predator individuals G to grow to the next size class i+1. According to the discretization of the McKendrick-von Foerster equation (see S1 File), we assume that G = e∑mp/(m − m), which means that the number of individuals transferred to size class i+1 is calculated by dividing the total biomass assimilated in size class i by the difference in individual masses between the neighboring size classes i and i+1. Then, where matrix ΔM+ = diag(m2 − m1,m3 – m2,…,m − m) describes the mass difference between the upper (i.e., subscript “+”) size class i+1 and size class i, and is the inverse of matrix ΔM+.

According to the growth process per unit time Δt in Fig 1A, G individuals from size class i–1 enter size class i, while G individuals in size class i leave size class i, resulting in a net change in the number of individuals in each size class i given of G – G. Additionally, there is an increase in the number of individuals in size class 1 corresponding to the recycling of carrion C from imaginary size class k+1, which is given by C/m1 (or Gm/m1). As the dimension of C is mass per unit time, division by m1 converts the mass into the number of resources. Thus, shifts in the number of individuals (Δn) caused by growth and the recycling of carrion from size class k+1 per unit time Δt are expressed as follows: where is a transfer matrix, with entries –1 on the diagonal and 1 on the subdiagonal, for transferring G to the next size class, and is a matrix to transfer carrion from size class k+1 to be recycled in size class 1.

In addition to the assimilated biomass in Eq (5), the total undigested (unassimilated) materials E per unit time Δt of each size class i is given by (1−e)∑mP. This can be equivalently written in matrix form as (I−E)NVSMn, where I is a unit matrix with ones on the diagonal and zeros elsewhere. For the recycling of these undigested materials into resources, a new transition matrix T1 to size class 1 is introduced: where Δn denotes the increase in the number of resources caused by recycling of undigested (unassimilated) materials per unit time Δt, and is the inverse of matrix M1 = diag(m1,m1,…,m1). Multiplying (I−E)NVSMn by converts E in Fig 1A, whose dimension is mass per unit time, into number of resources; multiplication by T1 then transfers all undigested materials from each size class i to size class 1.

In conclusion, by combining Eqs (4), (7) and (10), we obtain Δn as the shift in numbers obtained from the predator–prey interaction process per unit time Δt:

Reproduction

As depicted in Fig 1B, reproduction occurs only in the size range [m3,m]. We assume that the mean proportions of females and males in each size class i are α and 1−α, respectively. The reproduction rate F is given by , where c is the total biomass of eggs produced by unit weight and r is the scaling exponent, which does not differ between sexes.

Newborn individuals (eggs or larvae) are placed in size class 2 and excretive matter (semen) is recycled into size class 1 as newly input resources, thereby causing an increase in the numbers of individuals in size classes 1 and 2 as follows where A = diag(α1,α2,…,α) is a sex ratio matrix, is the reproduction rate matrix, T2 is a new matrix for transferring eggs into size class 2, where and is the inverse of matrix M2 = diag(m2,m2,…,m2).

Assuming that B individuals per unit time Δt are transferred from size class i to the next-lower size class i–1 owing to mass decrease by reproduction and is calculated using . The loss in the number of individuals in size class i as a result of transfer to size class i–1 and the gain in the number of individuals in size class i as a result of transfer from size class i+1 lead to a net change −B+B in the number of individuals in each size class i be given by where Δn denotes the shift in the number of individuals as a result of this process and is the inverse of matrix ΔM_ = diag(m1−m0,m2−m1,…,m−m) (where m0 = 0) that describes the mass difference between size class i and the lower (hence, the subscript “−”) size class i–1.

Therefore, the column vector Δn is given by adding Eqs (13) and (15):

Metabolism

Individuals must pay metabolic costs to maintain body structure and activity. Standard metabolic rate is recognized as an allometric function of body mass: vmΔt [24]. As shown in Fig 1C, the model of metabolism is similar to that of reproduction. Although resources do metabolize, the metabolites are eventually recycled back to their own size class, resulting in no variation in the number of individuals. Therefore, for simplicity only the metabolic costs of the consumers’ group are considered in our model. Recycling of these metabolites leads to an increase in the number of individuals in size class 1 as follows: where is the metabolic rate matrix whose diagonal elements are the metabolic rates W of the individuals in size class i.

If we assume that number of individuals transferred from size class i to the next-lower size class i–1 is given by , then the net change −L+L in the number of individuals in each size class i per unit time Δt is given by where Δn denotes the shift in the number of individuals as a result of this process.

Therefore, the resulting column vector Δn is the sum of Eqs (17) and (18):

Non-predation and fishing mortalities

The causes of non-predation mortality include disease, competition, senescence, and other natural factors. McCoy and Gillooly [19] stated that instantaneous natural mortality rate scales allometrically with body mass as μ0m−. However, the predation mortality for the largest individuals in our model is zero, while non-predation mortality is also naturally very low, resulting in an excessive accumulation of individuals. Increasing or decreasing the mortality of top-predators can cause a top-down cascade that leads to a change in the size distribution of their prey individuals and produces a wave-like biomass spectrum [26].

To obtain a linear size spectrum, we modulate the mortality of large individuals by introducing a density-dependent mortality as in Benoît and Rochet [4] but with the density-dependent mortality increasing with the ratio between the number of individuals and the carrying capacity. We assume that the carrying capacity is size-based and given by km−, where k is the carrying capacity factor indicating the maximum number of individuals that the community can support per unit mass and f is the scaling exponent. Consequently, the non-predation mortality of each individual in size class i becomes , where z0 = μ0/k and z = d−f.

Thus, the decreased number of individuals D owing to non-predation mortality (Δn) in each size class i per unit time Δt is where the product , whose diagonal element gives the non-predation mortality rate of each individual in size class i.

The recycling of carrion from each size class i yields new production (Δn) in size class 1: where the entries of the product NDMn give the biomass of the carrion of each size class i.

Fishing can affect community structures in marine ecosystems. Several community size spectrum models and trait-based size spectrum models have been used to investigate the effect of fishing [4, 11–13], and fishing mortality can be easily included in this matrix model to examine the effects of different fishing patterns on the size spectrum. Similar to the process of non-predation mortality, the decrease in the number of individuals as a result of fishing is given as where H = diag(h1,h2,…,h) is the fishing mortality rate matrix, whose diagonal element gives the fishing mortality rate of each individual in size class i. The recycling of these catches (Δn) through human consumption into size class 1, as depicted in Fig 1D, is then expressed as

Finally, combining Eqs (20), (21), (22) and (23), we obtain:

Unified size-structured matrix model

The overall matrix model for expressing size-structured community dynamics can be obtained by adding Eqs (12), (16), (19) and (24):

Parameter settings for numerical simulation

The egg size of marine fish is fairly constant among species and lies predominantly within a narrow range around 1 mg [27]. Thus, we set the mass of the first consumer size class as m2 = 0.001 g. To obtain a smooth line, we use a small value for the mass ratio of successive size classes, δ. The size range of imaginary resources [m−,m0] is also discretized using the ratio δ: m = mδ, where x = −s, −s + 1…,0. This range roughly extends from m− = 1.36×10−10 g to m0 = 4.44×10−4 g, which also covers very small phytoplankton cells [28]. To cover the size range of large adult fish in the ocean, the size axis is discretized over 50 size classes ranging from m1 = 6.67×10−4 g to m = 2.83×105 g. We assume that the total biomass of the community is 0.5 g/m-3 and run the model starting from an initial condition until it reaches a stationary state for each numerical simulation. The simulations are performed in Matlab (see code in S2 File). The steady-state conditions for size spectra are achieved within a 40-year model run (i.e., Eq (25)). It has been determined that the results do not depend on the initial number of individuals within each size class.

The other parameter values used in the basic model simulation case are listed in Table 1. Most of the parameter values are based on previously published size-based models. Parameters such as the exponents of the searching and reproduction rates are given as 3/4 based on the fact that metabolic rate scales to the 3/4 power of body size. The value of the reproduction rate factor c is set such that the ratio between the reproduction and assimilation rates falls within the reference range [6]. The values of factor k and the exponent f of the carrying capacity are derived from modified values of z0 and z in Benoît and Rochet [4].

Increasing the mass ratio δ of successive size classes leads to a decrease in the number of size classes; if the number is reduced by half, then it is necessary to increase the number of individuals in each size class to twice the number of individuals in the basic case to keep the total biomass in the community constant. In this case, the carrying capacity factor k should also be doubled to support the abundance of individuals. As shown in S1 Fig, varying the value of δ does not significantly change the size spectrum. Only the egg size class (i.e., size class 2) and its adjacent size classes depart slightly from the linear size spectrum because the total energy allocated to reproduction from consumers decreases with increasing value of δ. The size spectra obtained with different values of time step Δt are exactly the same (S1 Fig), implying that the time step does not change the dynamics of size spectra. A mass ratio δ of successive size classes larger than 1.5 will allow the model to run at a daily or a larger time step (S1 Fig).

Size-dependency of rates in each process

Using the parameter values in Table 1 results in a mass-specific growth and reproduction rates that decrease with body mass with slopes close to –1/4 (S2 Fig) on a log–log scale. The log–log plot of predation mortality rate versus body size is nearly linear, with a slope close to –1/4 within the size range between 10−3 and 103 g. For individuals larger than 103 g, predation mortality decreases sharply with body size owing to the decline in predator individuals (S2 Fig). The logs of non-predation mortality rate and body size show a nearly linear relationship, with higher and lower mortality rates for smaller and larger individuals, respectively, although with a slope that is slightly shallower than –1/4 (S2 Fig). This may be caused by carrying capacity effects.

Results

Steady state

Using the reference values of the parameters in Table 1 results in a linear size spectrum with a slope of approximately –1, which is consistent with values reported in empirical studies [4, 6] (Fig 2).

Fig 2

Linear size spectrum.

Simulated size spectrum at steady state (dotted line) using reference parameters in Table 1 and a regression line as a function of log body size (solid line; log10[abundance] = −0.9959×log10[body mass] − 2.006).

Examples of numerical simulation

In the application shown in Fig 2, we simulated the size spectrum with the constant preferred PPMR β given in Table 1. However, field data show that there is a great variability in β with predator mass [29, 30]. Therefore, the effect of varied β () was investigated for τ ≠ 0. For a plausible range of τ (−0.2 ≤ τ ≤ 0.2), the size spectra are essentially identical to the basic case spectrum (Fig 3A). An extreme range of τ (−0.9≤τ≤0.9), which apparently exceeds the plausible range obtainable in an actual situation, was also examined to obtain a clear contrast in terms of the effect of τ. For υ = 105 and τ = −0.9, for which β ranges from approximately 5×107 for small predators of mass 10−3 g to approximately 1.0 for the largest predators of mass 3×105 g, static waves appear along the size spectrum (Fig 3B). By contrast, for υ = 105 and τ = 0.9, for which β increases with predator body mass and ranges from approximately 2×102 for small predators of 10−3 g to approximately 8×109 for the largest predators, the resulting size spectrum shows a slightly steeper slope (Fig 3B). Under various combinations of υ and τ, the wavelengths of the size spectra vary (S3 Fig). For all cases the slopes range from approximately –0.98 to –1.16 but are close to that of the basic case shown in Fig 2 (ca. –1). Thus, except for the wavelengths along the size spectrum, PPMR does not significantly affect the size spectrum.

Fig 3

Simulated size spectra at equilibrium with different parameter values.

Simulated size spectra using (a) plausible combinations of the factor υ and exponent τ of PPMR, (b) extreme combinations of υ and τ, different values of (c) the reproduction rate factor c, (d) reproduction rate exponent r, (e) metabolic rate exponent q with a small metabolic rate factor v, and (f) exponent q with a large factor v. The vertical dotted lines mark wave troughs or peaks and the horizontal distance between two adjacent vertical dotted lines is the wavelength. The horizontal arrows indicate the length of 2log10β.

Reduction in the value of c or r in the reproductive rate leads to a decrease in reproductive allocation and an eventual accumulation of large individuals, which in turn can cause the population of their most preferred prey to decrease through predation as a top-down effect, resulting in the alternate propagation of opposed effects to smaller size ranges. As a result, waves propagate through the size spectrum (S4 Fig). Size spectra with static waves emerge when the system reaches equilibrium conditions (Fig 3C and 3D). The wavelength depends on PPMR [31], with observed wavelengths approximately equal to or less than 2log10β on a log–log scale (Fig 3C and 3D). Conversely, increasing the value of c or r leads to a decrease in the number of large individuals and a cut-off of the size spectrum (Fig 3C and 3D), which arises because large individuals convert more energy into reproduction, leaving less energy available for growth.

The standard metabolic rate vm is smaller than both the predation and reproduction rates and therefore might have a weak effect on the size spectrum. To show how the size spectrum can be affected by the standard metabolic rate, we set the factor v at different levels for comparison. For low values of v (Fig 3E), the standard metabolic rate of each individual in size class i is very small, and therefore changing the exponent q has no apparent effect on the size spectrum. For large values of v, however, varying q changes the shape of the size spectrum (Fig 3F). A size spectrum with static waves with wavelengths less than 2log10β appears when q decreases from the reference value 0.75 to a value of 0.25. Before reaching an equilibrium state, waves travel from small to large individuals (S4 Fig). An increase in q produces a size spectrum with a cut-off in the range of large body masses.

Effect of fishing on the size spectrum

Two types of fishing patterns were compared: (i) an unbalanced fishing pattern in which individuals larger than 104 g are caught with a constant fishing mortality q0F0, where q0 is the catchability coefficient and F0 is the fishing intensity coefficient; and (ii) a balanced fishing pattern in which individuals larger than 10 g are caught in proportion to their natural productivity determined by the metabolic size-scaling rules [18, 32], i.e., the fishing mortality rate for individuals in this size range is given by . The fishing mortality rate matrices corresponding to fishing patterns (i) and (ii) are given by H1 = diag(0,…,0, q0F0,…,q0F0) and , respectively (S5 Fig).

The unbalanced fishing pattern results in a more severe cut-off when the fishing intensity coefficient F0 increases from 5 yr-1 to 50 yr-1; by contrast, the balanced fishing pattern does not result in a cut-off of the spectrum even when the fishing intensity coefficient increases (Fig 4). Yields from the balanced fishing pattern are much higher than those of the unbalanced fishing pattern (S6 Fig).

Fig 4

Predicted impacts of fishing.

Effects of unbalanced and balanced fishing on the size spectrum at different fishing intensity coefficient levels: (a) F0 = 5 yr−1, (b) F0 = 10 yr−1, (c) F0 = 30 yr−1, and (d) F0 = 50 yr−1.

Discussion

We developed a size-structured matrix model to express community dynamics in which each process is modeled by simple body size power functions. The model can be seen as a discretization of the Mckendrick-von Foerster partial differential equation that produces a diffusion of growth, i.e., during each time step, several individuals grow up to the next-higher size class while other individuals remain in the same size class. Thus, the basic difference between our matrix model and the continuous Mckendrick-von Foerster equation lies in its ability to incorporate growth variation in individuals within a size class; in a real process, this diffusion is appropriate because individuals with the same size in a given population often differ in growth rate [33]. The main advantage of our model lies in its flexibility in modeling each process. In addition, the model can be run on a daily or larger time step Δt by increasing the mass ratio δ of successive size classes. In this case, the running time for numerical simulations required by our matrix model will be much shorter than that of PDE models as a result of the reduced number of size classes and loop iterations.

Boundary conditions imposed in upper and lower size classes can affect the shapes of size spectra. For the upper boundary, the introduction of senescence mortality in previous models [14, 34–36] led to a sharp increase in the mortality of large individuals. As a result, large size classes depart from the linear size spectrum, resulting in a downwardly curved stationary size spectrum and the emergence of a cut-off. However, field data support the assertion that natural mortality for organisms is negatively correlated with body mass [19, 37]. We therefore adopt an allometric natural mortality rate against body mass [11] based on empirical information. Our model is developed with a type I functional response, which can produce an unlimited growth of predators and an unrealistic increase in the number of these individuals. Introducing a density-dependent effect to the allometric natural mortality rate can help control the increase in predator individuals, thereby improving the stability of size spectra.

For the lower boundary, the growth and reproduction of background resources based on the processes of nutrient uptake and photosynthesis are not explicitly modeled in our study; nevertheless, our model can successfully reproduce linear size spectra with slopes close to –1 (Fig 2). The size spectrum produced by the basic case does not depart from a linear shape for nearly all size classes, possibly because we introduced an imaginary size range to allow consumers to have a common food availability level.

Under the mass balance assumption, variation in the biomass (and abundance) of background resources is controlled by the difference between the rates of recycling of biomass egested from consumers and consumption of resource biomass by predators. In our model, the community can reach equilibrium conditions even under extreme parameter values, although size spectra with cut-offs and/or static waves are frequently observed in such cases (Fig 3). If the biomass of resources is assumed to be constant [14, 31], infinite oscillations occur when the non-predation mortality rate is set to zero, while only transient oscillations emerge in this case under the mass balance assumption (S7 Fig). Under the constant resource assumption, there are some cases in which the total biomass of the community will continue to decrease, eventually leading to a collapse of the community under some extreme parameter value assumptions, e.g., r = 0.25 (S8 Fig). Owing to the mass balance assumption, the effect of extreme parameter values on the size spectrum can be alleviated. If model parameters are varied over reasonable ranges, similar results can still be obtained (S8 Fig). However, because carbon is exchanged and recycled throughout the biosphere and the mass balance in a local community might not always be retained, it will be important to examine the consequences of various types of assumptions on background resources in future studies that explore the dynamics of real marine communities.

Maury et al. [6] found that their model produces unstable oscillations when the reproductive investment of females is extremely low. Increasing the value β of the most preferred PPMR or decreasing the feeding kernel width σ can also make the size spectrum unstable and generate traveling waves (e.g., [31]). Unlike these previous studies, neither infinite oscillations nor traveling waves appeared in our study. In addition to the mass balance assumption discussed above, this may be because of the presence of growth diffusion in our model, which can help stabilize the size spectrum [14].

Thiebaux and Dickie [38] mathematically modeled a predator biomass spectrum with a parabolic dome and predicted that the related prey biomass spectrum would also become a parabola as a result of energy transfer between trophic levels, which explained their corresponding wave-like size spectra. Our results suggest that such waves are caused by the accumulation of large individuals, as a top-down cascade can be stimulated by an increase in the number of predators with mass m in the largest size class suppressing their prey individuals with mass m/β, which in turn releases individuals with mass from predation. This should result in a wavelength approximately equal to 2log10β on a log–log scale; however, our results primarily produced wavelengths less than 2log10β. This may have been caused by growth rates that differed between predators and their preferred prey: in this case, because small individuals travel faster than large individuals along the size spectrum, the distance between the peak of one wave and the next-lower trough would become smaller than PPMR.

Wave-like size spectra have also been observed in nature. For example, Gaedke [39] reported wave-like oscillations in the seasonal distribution of the biomass of individuals ranging from bacteria to crustaceans. These results imply that a simple linear size spectrum will not be sufficient to describe the dynamics of actual marine communities. Wave-like oscillations also need to be recognized as important indicators of ecosystems, and if simple linear regression analysis alone is relied upon, the presence of cut-offs and/or waves might result in inadequate slope estimates. The shape of a size spectrum should therefore be looked at closely in revealing the impacts of human activities and environmental changes.

To show the applicability of our model, a short example concerning the fishing effects on a marine community was enumerated (see the last part of Results). The present matrix model can be easily expanded to enable investigation of the effects of various other situations, e.g., stochastic or periodic environmental changes such as water temperature and productivity, which vary with season [40]. Other processes such as seasonal reproduction and long-term environmental changes can also be coupled to the model to investigate the response of marine community size spectra at large scales.

The model can be easily modified to take species identity into consideration and include population demography. To achieve this, we can replace the column vector n(t) with a new column vector, e.g., n((t) to represent the structure of species 1. For instance, in the predator–prey interaction process the total number of individuals of species 1 in each size class eaten by predators of other species h () per unit time Δt can be expressed as , where N( is the abundance matrix of species h and θ( is a matrix representing the interaction strength between species h and 1. The elements in the interaction strength matrix can be freely varied for application to several cases, e.g., in which species’ prey preferences vary by prey species and size. Likewise, species-specific reproduction can be modeled by setting the maturation and maximum sizes in each species’ abundance vector (e.g., n((t) = (n((t), n((t),…, n((t),…, n((t))T). The reproduction rate for each species can then be given by AR(n(, where R( is a reproduction rate matrix of species h whose elements are modified to place the maturation, maximum, and offspring sizes within the proper size classes for the species. Other investigations on, for instance, the impacts of various fishing patterns for different species on marine communities can also be conducted through similar structural modifications of the present model. Consequently, a three-dimensional matrix model formed by three coordinate axes, “body mass” (x), “species” (y), and “abundance” (z), is envisioned.

An important goal for future work is to apply the model to real ecosystem. To produce reasonable results, the model and parameter values in the present study can be calibrated for different ecosystem types based on recorded data related to the vital rates, life histories, and behaviors of individuals within a specific community. Determining appropriate lower and upper boundary conditions will be especially important in calibrating the model to reproduce real-world dynamics. Real size distribution of communities obtained from research vessel data and other sources should be compared with the simulation results derived from the calibrated model through the use of sets of sufficiently powerful statistical measures.

Conclusion

In this paper, we demonstrated the construction of a community matrix model with several body size power functions from individual-level processes. The model can produce stable size spectra under various conditions. Growth diffusion, which plays a role in stabilizing size spectra, is allowed in our model. When parameter values change, the model produces static waves at equilibrium states as a result of top-down cascades or cut-offs at large size ranges caused by reduction in energy allocated to growth or high fishing mortality.

The matrix algorithms developed in this study can reduce the computational time needed to calculate size spectra. Increasing the mass ratio of successive size classes enables a daily or larger time step without affecting the results. As a result of its flexibility, the model can be easily extended or modified to explore community- or population-level dynamics.

Influence of mass ratio and time step on the size spectrum at steady state.

Normalized size spectra, in which the abundance is divided by the logarithmic width of each size class, at steady states for different values of (a) mass ratio δ of successive size classes and (b) time step Δt.

(TIF)

Click here for additional data file.

Size-dependent relationships used to parameterize the model.

(a) Decrease in mass-specific growth, reproduction, and metabolic rates with increasing size. (b) Decrease in predation mortality, and non-predation mortality rates with increasing size.

(TIF)

Click here for additional data file.

Effects of varying PPMR on the size spectrum.

Simulated size spectra at steady state for different combinations of factor υ and exponent τ of PPMR.

(TIF)

Click here for additional data file.

Examples of transient oscillations prior to reaching equilibrium.

When (a) the reproduction rate factor c is set to 0.1 and (b) the metabolic rate factor v and exponent q are set to 5 and 0.25, respectively, the model produces transient oscillations until the system reaches a stationary state. Dotted lines represent the waves propagating through the spectrum, while solid lines represent the stationary states (i.e., static waves).

(TIF)

Click here for additional data file.

Fishing mortality.

Relationship between body mass and fishing mortality for four different fishing intensity coefficients F0 under unbalanced and balanced fishing patterns.

(TIF)

Click here for additional data file.

Yields from the two fishing patterns under different fishing intensity coefficients.

(TIF)

Click here for additional data file.

Comparison of equilibrium and non-equilibrium states.

When the non-predation mortality rate factor μ0 is set to zero, (a) static waves appear under mass balance consumption, while (b) traveling waves emerge under constant resource assumption. Solid line represents the steady state and dotted lines represent waves that change over time.

(TIF)

Click here for additional data file.

Examples of simulated size spectra with a constant background resource.

Different values of (a) the reproduction rate factor c, (b) the reproduction rate exponent r, (c) the metabolic rate exponent q with a small metabolic rate factor v, and (d) the exponent q with a large factor v are used for simulations. Solid lines represent size spectra obtained from running the system to steady state, while dotted lines represent collapsed communities that will never reach equilibrium.

(TIF)

Click here for additional data file.

Calculating the number of individuals G and B.

(DOCX)

Click here for additional data file.

Matlab code.

(DOCX)

Click here for additional data file.