A numerical framework for computing steady states of structured population models and their stability

  • Received: 01 January 2016 Accepted: 01 January 2016 Published: 01 August 2017
  • MSC : Primary: 35B40, 92B05; Secondary: 65N40

  • Structured population models are a class of general evolution equations which are widely used in the study of biological systems. Many theoretical methods are available for establishing existence and stability of steady states of general evolution equations. However, except for very special cases, finding an analytical form of stationary solutions for evolution equations is a challenging task. In the present paper, we develop a numerical framework for computing approximations to stationary solutions of general evolution equations, which can \emph{also} be used to produce approximate existence and stability regions for steady states. In particular, we use the Trotter-Kato Theorem to approximate the infinitesimal generator of an evolution equation on a finite dimensional space, which in turn reduces the evolution equation into a system of ordinary differential equations. Consequently, we approximate and study the asymptotic behavior of stationary solutions. We illustrate the convergence of our numerical framework by applying it to a linear Sinko-Streifer structured population model for which the exact form of the steady state is known. To further illustrate the utility of our approach, we apply our framework to nonlinear population balance equation, which is an extension of well-known Smoluchowski coagulation-fragmentation model to biological populations. We also demonstrate that our numerical framework can be used to gain insight about the theoretical stability of the stationary solutions of the evolution equations. Furthermore, the open source Python program that we have developed for our numerical simulations is freely available from our GitHub repository (github.com/MathBioCU).

    Citation: Inom Mirzaev, David M. Bortz. A numerical framework for computing steady states of structured population models and their stability[J]. Mathematical Biosciences and Engineering, 2017, 14(4): 933-952. doi: 10.3934/mbe.2017049

    Related Papers:

    [1] Jordi Ripoll, Jordi Font . Numerical approach to an age-structured Lotka-Volterra model. Mathematical Biosciences and Engineering, 2023, 20(9): 15603-15622. doi: 10.3934/mbe.2023696
    [2] Weicheng Chen, Zhanping Wang . Stability of a class of nonlinear hierarchical size-structured population model. Mathematical Biosciences and Engineering, 2022, 19(10): 10143-10159. doi: 10.3934/mbe.2022475
    [3] Marcin Choiński, Mariusz Bodzioch, Urszula Foryś . A non-standard discretized SIS model of epidemics. Mathematical Biosciences and Engineering, 2022, 19(1): 115-133. doi: 10.3934/mbe.2022006
    [4] Thanh Nam Nguyen, Jean Clairambault, Thierry Jaffredo, Benoît Perthame, Delphine Salort . Adaptive dynamics of hematopoietic stem cells and their supporting stroma: a model and mathematical analysis. Mathematical Biosciences and Engineering, 2019, 16(5): 4818-4845. doi: 10.3934/mbe.2019243
    [5] Cristeta U. Jamilla, Renier G. Mendoza, Victoria May P. Mendoza . Explicit solution of a Lotka-Sharpe-McKendrick system involving neutral delay differential equations using the r-Lambert W function. Mathematical Biosciences and Engineering, 2020, 17(5): 5686-5708. doi: 10.3934/mbe.2020306
    [6] Elvira Barbera, Giancarlo Consolo, Giovanna Valenti . A two or three compartments hyperbolic reaction-diffusion model for the aquatic food chain. Mathematical Biosciences and Engineering, 2015, 12(3): 451-472. doi: 10.3934/mbe.2015.12.451
    [7] Karl Peter Hadeler . Structured populations with diffusion in state space. Mathematical Biosciences and Engineering, 2010, 7(1): 37-49. doi: 10.3934/mbe.2010.7.37
    [8] Yansong Pei, Bing Liu, Haokun Qi . Extinction and stationary distribution of stochastic predator-prey model with group defense behavior. Mathematical Biosciences and Engineering, 2022, 19(12): 13062-13078. doi: 10.3934/mbe.2022610
    [9] Guangrui Li, Ming Mei, Yau Shu Wong . Nonlinear stability of traveling wavefronts in an age-structured reaction-diffusion population model. Mathematical Biosciences and Engineering, 2008, 5(1): 85-100. doi: 10.3934/mbe.2008.5.85
    [10] Jessica L. Hite, André M. de Roos . Pathogens stabilize or destabilize depending on host stage structure. Mathematical Biosciences and Engineering, 2023, 20(12): 20378-20404. doi: 10.3934/mbe.2023901
  • Structured population models are a class of general evolution equations which are widely used in the study of biological systems. Many theoretical methods are available for establishing existence and stability of steady states of general evolution equations. However, except for very special cases, finding an analytical form of stationary solutions for evolution equations is a challenging task. In the present paper, we develop a numerical framework for computing approximations to stationary solutions of general evolution equations, which can \emph{also} be used to produce approximate existence and stability regions for steady states. In particular, we use the Trotter-Kato Theorem to approximate the infinitesimal generator of an evolution equation on a finite dimensional space, which in turn reduces the evolution equation into a system of ordinary differential equations. Consequently, we approximate and study the asymptotic behavior of stationary solutions. We illustrate the convergence of our numerical framework by applying it to a linear Sinko-Streifer structured population model for which the exact form of the steady state is known. To further illustrate the utility of our approach, we apply our framework to nonlinear population balance equation, which is an extension of well-known Smoluchowski coagulation-fragmentation model to biological populations. We also demonstrate that our numerical framework can be used to gain insight about the theoretical stability of the stationary solutions of the evolution equations. Furthermore, the open source Python program that we have developed for our numerical simulations is freely available from our GitHub repository (github.com/MathBioCU).


    1. Introduction

    Many natural phenomena can be formulated as the differential law of the development (evolution) in time of a physical system. The resulting differential equation combined with boundary conditions affecting the system are called evolution equations. Evolution equations are a popular framework for studying the dynamics of biological populations. For example, they have proven useful in understanding the dynamics of biological invasions [51], bacterial flocculation in activated sludge tanks [7], and the spread of parasites and diseases [30]. Since many biological populations converge to a time-independent state, many researchers have used theoretical tools to investigate long-term behavior of these models. Analytical and fixed point methods have been used to show the existence of stationary solutions to size-structured population models [44,29] and semigroup theoretic methods have been used to investigate the stability of these stationary solutions [28,43,4]. For many models in the literature, the principle of linearized stability [55,34] can be used to show that the spectral properties of the infinitesimal generator (IG) of the linearized semigroup determines the stability or instability of a stationary solution. Moreover, using compactness arguments, spectral properties of the infinitesimal generator can be determined from the point spectrum of the IG, which in turn can be written as the roots of a characteristic equation.

    Despite this theoretical development, the derived existence and stability conditions are oftentimes rather complex, and accordingly the biological interpretation of these conditions can be challenging. To overcome this difficulty several numerical methods for stability analysis of structured population models have been developed [12,27,37,21]. For instance, Diekmann et al. [37,22,21] present a numerical method for equilibrium and stability analysis of physiologically structured population models (PSPM) or life history models, where individuals are characterized by a (finite) set of physiological characteristics (traits such as age, size, sex, energy reserves). In this method a PSPM is first written as a system of integral equations coupled with each other via interaction (or feedback) variables. Consequently, equilibria and stability boundary of the resulting integral equations are numerically approximated using curve tracing methods. Later, de Roos [19] presented the modification of the curve tracing approach to compute the demographic characteristics (such as population growth rate, reproductive value, etc) of a linear PSPM. For a fast and accurate software for theoretical analysis of PSPMs we refer interested reader to a software package by de Roos [20]. An alternative method for stability analysis of physiologically structured population models, developed by Breda and coworkers [11,13,10], uses a pseudospectral approach to compute eigenvalues of a discretized infinitesimal generator. This method (also known as infinitesimal generator (IG) approach) has been employed to produce bifurcation diagrams and stability regions of many different linear evolution equations arising in structured population modeling [13,14,15]. Unfortunately, not all structured population models fit into the framework of PSPMs and thus there is a need for a more general numerical framework for asymptotic analysis of structured population models.

    In this paper we develop a numerical framework to guide theoretical analysis of structured population models. We demonstrate that our methodology can be used for numerical computation and stability analysis of positive stationary solutions of both linear and nonlinear size-structured population models. Moreover, we illustrate the utility of our framework to produce approximate existence and stability regions for steady states of size-structured population models. We also provide an open source Python program used for the numerical simulations in our GitHub repository [42]. Although, the examples presented in this paper are size-structured population models, in Section 2, we show that the framework is applicable to more general evolution equations.

    The main idea behind the numerical framework is first to write a structured population model in the form of an evolution equation and then use the well-known Trotter-Kato Theorem [53,35] to approximate the infinitesimal generator of the evolution equation on a finite dimensional space. This in turn allows one to approximate solutions (or spectrum) of the evolution equation with the solution (or spectrum) of system of differential equations. Consequently, we approximate the stationary solutions of an actual model with stationary solutions of the approximate infinitesimal generator on a finite dimensional space. Approximate local stability of the approximate steady states are then computed from the spectrum of the Jacobian of ODE system evaluated at this steady states. Our method is similar to the IG approach in [13,14,15], in a sense that we also approximate infinitesimal generator and analyze the spectrum of the resulting operator to produce existence and stability regions. However, in contrast to IG approach, our framework also computes actual steady states and is easily applied to nonlinear evolution equations arising in structured population dynamics.

    The rest of the paper is structured as follows. We describe the theoretical development of our framework for general evolution equations in Section 2.1. Note that readers with more biological background can skip Section 2.1 and directly jump into the application of the framework in Section 2.3. In Section 2.3, we illustrate the convergence of the approximation method by applying it to linear Sinko-Streifer model, for which the exact form of the stationary solutions is known. To further illustrate the utility of our approach, in Section 3, we apply our framework to a nonlinear size-structured population model (also known as population balance equations in the engineering literature) described in [5,9]. Moreover, in Section 2.2, we show that approximate local stability conditions for a stationary solution can be derived from the spectral properties of the approximate infinitesimal generator. Finally, we conclude with some remarks and a summary of our work in Section 5.


    2. Numerical framework

    In this section, we demonstrate our numerical methodology for general evolution equations. We first present the numerical scheme used to approximate the infinitesimal generator of a semigroup. Subsequently, in Section 2.3, we illustrate the convergence of our approach by applying it to linear Sinko-Streifer equations, for which exact stationary solutions are known.


    2.1. Infinitesimal generator approximation

    Given a Banach space X, consider an abstract evolution equation,

    ut=F(u),u(0,)=u0X, (1)

    where F:D(F)XX is some operator defined on its domain D(F) and u0 is an initial condition at time t=0. Note that any boundary condition belonging to a given partial differential equation can be included in the domain D(F). The solution to (1) can be related to the initial state u0 by some transition operator T(t) so that

    u(t,x)=T(t)u0(x).

    The transition operator T(t) is said to be a strongly continuous semigroup (or simply C0 -semigroup) if satisfies the following three conditions:

    1. T(s)T(t)=T(s+t) for all s,t0

    2. T(0)=I, the identity operator on X

    3. For each fixed u0X,

    limt0+T(t)u0u0=0.

    Moreover, showing that the operator F generates a C0-semigroup is equivalent to establishing well-posedness of the abstract evolution equation given in (1). In general, finding the explicit form of the transition operator is challenging. Hence, approximation methods, e.g. Yosida approximants, are used to study a more complicated evolution equation and the semigroups they generate. One of the famous approximation theorems is due to Trotter [53] and Kato [35] (see [36] for the classical literature on the approximation of generators of semigroups). The goal is to construct approximating generators Fn on the approximate spaces Xn such that C0-semigroups Tn() generated by Fn approximate the C0-semigroup T(t) generated by F. Although there are multiple ways to approximate the infinitesimal generator F, for our purposes we use the approximation scheme based on Galerkin-type projection of the state space X [6,32,2]. For the convenience of readers, we will summarize the approximation scheme here.

    Let Xn, n=1,2, be a sequence of subspaces of X with dim(Xn)=n and define projections πn:XXn and canonical injections ιn:XnX. Assume that the projections πn are bounded, i.e., there exists ˜M>0 such that

    πn˜M (A1)

    for all n=1,,n. Moreover, assume that

    limnπnv=v (A2)

    for all vX. Consequently, for each subspace Xn we choose basis elements {βni}ni=1 such that each element v of subspace Xn can be written in the form v=ni=1αiβni. Moreover, for each subspace Xn we define the bijective mappings pn:XnRn by

    pnv=(α1,,αn)T

    and the norm on Rn by

    vRn=p1nvX.

    Consequently, we define bounded linear operators Pn:XRn and En:RnX by

    Pnv=pnπnv,vX (2)

    and

    Enz=ιnp1nz,zRn, (3)

    respectively. Finally, we define approximate operators Fn on Rn by

    Fn(z)=PnF(Enz) (4)

    for all zRn.

    Accordingly, the Trotter-Kato Theorem states that the semigroup generated by the discrete operator Fn converges to the semigroup generated by the infinitesimal generator F. For the convenience of the reader, we state the theorem here and refer readers to [32] for a proof.

    Theorem. (Trotter-Kato) Assume that (A1) and (A2) are satisfied. Let (T(t))t0 and (Tn(t))t0, nN, be strongly continuous semigroups on X and Rn with generators F and Fn, respectively. Furthermore, assume that they satisfy the estimate

    T(t)X,Tn(t)RnMewtfor all t0,nN,

    for some constants M1,wR. Then the following are equivalent:

    1. There exists a λ0ρ(F)ni=1ρ(Fi) such that for all vX

    En(λ0InFn)1Pnv(λ0IF)1vX0as n.

    2. For all vX and t0,

    EnTn(t)PnvT(t)vX0

    as n, uniformly on compact t intervals.

    In general, one establishes the first statement for a Trotter-Kato approximation and then uses the second statement to approximate an abstract evolution equation on a finite dimensional space. In their paper, Ito and Kappel [32] present the standard ways to establish the first statement of the theorem (see also [6,2,1]). Therefore, here we assume that for a particular problem the first statement in the theorem has already been established and thus the evolution equation in (1) can be approximated by the following system of ODEs,

    un(t)=Fn(un(t)),un(0)=Pnu(0,). (5)

    Consequently, the solution of the IVP is mapped onto the infinite dimensional Banach space X and one has the following convergence

    limnEnun(t)uX=0 (6)

    for t in compact intervals.

    In general, finding explicit stationary solutions of abstract evolution equations is a challenging task. Conversely, many efficient root finding methods have been developed for finding steady states of a system of ODEs. For large-scale nonlinear systems, many efficient methods have been developed as well. Hence, we propose a numerical framework that utilizes those efficient root finding methods to approximate steady state solutions of general evolution equations. The idea is to use an efficient and accurate root finding method to approximate a stationary solution of the evolution equation (1) with the stationary solutions of the IVP in (5). Thus, as a consequence of the Trotter-Kato Theorem, the steady states of (5) converge to the steady states of (1) as n.


    2.2. Stability of stationary solutions

    Studying the asymptotic behavior of solutions is a fundamental tool for exploring the evolution equations which arise in the mathematical modeling of real world phenomena. To this end, many mathematical methods have been developed to describe long-term behavior of evolution equations. For instance, for long-time behavior of linear evolution equations, linear semigroup theoretic methods can be used to formulate physically interpretable conditions. Furthermore, for nonlinear evolution equations, the principle of linearized stability can be used to relate the spectrum of the linearized infinitesimal generator to the local stability or instability of the stationary solution. Nevertheless, investigating the spectrum of the linearized infinitesimal generator is cumbersome and requires advanced functional analysis techniques. In contrast to general evolution equations, the asymptotic behavior of ordinary differential equations are determined by the eigenvalues of the Jacobian and well-studied. Hence, in this section we demonstrate that the approximation scheme, presented in Section 2.1, can also be used to give some insights about the stability of stationary solutions of the general evolution equations.

    Stability results discussed in this section are not in a traditional Lyapunov sense. In particular, since stationary solutions discussed in this paper are only approximations to actual stationary solutions, the stability results only hold for finite time intervals. Therefore, we refer to this kind of stability as approximate local stability of stationary solutions as this stability is deduced from numerical approximation of an evolution equation. In mathematical terms, the local numerical stability of a stationary solution is defined as follows.

    Definition 2.1. Stationary solution u of an abstract evolution (1) is called approximately locally stable, if for every closed finite time interval [0,tf] and ϵ>0, there exists δ>0 such that a unique solution of (1), u(t,x), with initial condition u0 fulfilling u0uX<δ satisfies

    u(t,)uX<ϵ (7)

    for all t[0,tf].

    Having the required definition in hand, we now prove the following stability result for general evolution equations.

    Corollary 2.2. Let u denote a stationary solution of the abstract evolution (1) and JA(un) denote the Jacobian of the approximate system of ODEs defined in (5). If for all sufficiently large n the eigenvalues of JA(Pnu) are strictly negative, then u is approximately locally stable in the sense of Definition 2.1.

    Proof. Since the infinitesimal generator approximation scheme, presented in Section 2.1, is convergent, for every given ϵ>0 and finite time interval [0,tf] there exist nϵN such that for nnϵ,

    u(t,)Enun(t)X<ϵ/2 (8)

    for all t[0,tf] (where the bounded linear function En is defined as in (3)). Moreover, the eigenvalues of JA(Pnu) are strictly negative for all sufficiently large n. This in turn implies that PMu is a locally asymptotically stable solution of (5) for some Mnϵ. That is, for given ϵ>0 there is δ>0 such that

    uM(t,)PMuRM=EMuM(t,)uX<ϵ/2 (9)

    for all t>0 and for all u0 such that PMu0PMu)RM=u0uX<δ (see for instance [323]). Consequently, combining (8) and (9) yields

    u(t,)uXu(t,)EMuM(t,)X+EMuM(t,)uX<ϵ

    for all t[0,tf] and for all u0 such that u0uX<δ.

    We note that although the stability result of Corollary 2.2 holds for arbitrarily large finite time intervals, the Corollary does not guarantee Lyapunov stability of stationary solutions.


    2.3. Numerical convergence results

    To verify convergence of the proposed approximation scheme, we apply the framework to the linear Sinko-Streifer model [52] for which an exact form of the stationary solution is available. The model describes the dynamics of single species populations and takes into account the physiological characteristics of animals of different sizes (and/or ages). The mathematical model reads as

    ut=G(u)=(gu)xμu,t0,0x¯x< (10)

    with a McKendrick-von Foerster type renewal boundary condition at x=0

    g(0)u(t,0)=¯x0q(y)u(t,y)dy

    and initial condition

    u(0,x)=u0(x).

    The variable u(t,x) denotes the population density at time t with size class x. The population is assumed to have a minimum and a maximum size 0 and ¯x<, respectively. The function g(x) represents the average growth rate of the size class x and the coefficient μ() represents a size-dependent removal rate due to death or predation. The renewal function q() represents the number of new individuals entering the population due to birth.

    Setting the right side of the equation (10) to zero and integrating over the size on (0,x) yields the exact stationary solution

    u(x)=1g(x)exp(x0μ(s)g(s)ds)¯x0q(y)u(y)dy. (11)

    Multiplying both sides of (11) by q(x) integrating over the size on (0,¯x), we obtain a necessary condition for existence of a stationary solution,

    1=¯x0q(x)g(x)exp(x0μ(s)g(s)ds)dx. (12)

    The convergence of the approximation scheme presented in Section 2.1 for Sinko-Streifer models has already been established in [6]. Using the basis functions for n-dimensional subspace Xn of the state space X=L1(0,¯x) are defined as

    βni(x)={1;xni1<xxni;i=1,,n0;otherwise

    for positive integer n with {xni}ni=0 a uniform partition of [0,¯x], and Δx=xnjxnj1 for all j. The functions βn form an orthogonal basis for the approximate solution space

    Xn={hX|h=ni=1αiβni,αiR},

    and accordingly, we define the orthogonal projections πn:XXn

    πnh(x)=nj=1αjβnj(x),whereαj=1Δxxnjxnj1h(x)dx.

    Moreover, since the evolution equation defined in (10) is a linear partial differential equation, the approximate operator Gn on Rn is given by the following n×n matrix

    [Gn]ij={1Δxg(xni)μ(xni)+q(xnj)i=j=1q(xnj)i=1,j21Δxg(xni)μ(xni)i=j21Δxg(xni1)i=j+1,j10otherwise. (13)

    At this point, one can use numerical techniques to calculate zeros of the linear system

    Gnun=0. (14)

    For the purpose of illustration, we set the model rates to

    q(x)=a(x+1),g(x)=b(x+1),μ(x)=c. (15)

    Plugging this rates into the necessary condition (12) yields

    a={ln2bb=c(bc)/(21c/b1)bc.

    For b=c this is a curve in 3D coordinate system and when bc the surface is illustrated in Figure 1a. To illustrate the utility of our approach, we used the numerical scheme described in this section to generate existence and stability regions for Sinko-Streifer model. Particularly, the interval (a,b,c)[0,1]×[0,1]×[0,1] is discretized with Δa=Δb=Δc=0.01. Consequently, we checked for the existence of a positive steady state at each of these discrete points. Since the approximate operator Gn is an n×n matrix, we can check if Gn is a singular matrix using standard tools. The resulting existence region is depicted in Figure 1a forming a nontrivial three-dimensional curve (red) on the surface defined by the necessary condition for existence of the steady states. Moreover, the existence and stability regions of the Sinko-Streifer model coincide for the chosen model rates in (15).

    Figure 1. Results-of-the Results of the numerical simulations. a) a,b,c values satisfying the necessary condition (12), form a 3D surface (blue surface). Steady states of the Sinko-Streifer model only exist on the red line b) Comparison of exact stationary solution (for the point marked with red star in Figure 1a) with approximate stationary solution for n=100. c) Absolute error between exact stationary solution and approximate stationary solution decays linearly as the dimension of approximate subspaces Xn increase.

    For the purpose of illustration, we arbitrarily choose b=0.5 and a=c=1 (marked with a red star in Figure 1a) as for these values a positive steady state of the Sinko-Streifer model exists. Figure 1b indicates that even for n=100 the fit between approximate and actual stationary solution is satisfactory (the infinity norm of the error is 0.004). Moreover, Figure 1c illustrates that as the dimension of the approximate space Xn increases the absolute error between the exact stationary solution and the approximate stationary solution decreases. Furthermore, the numerical algorithm has a linear convergence rate. This is due to the fact that we chose zeroth order functions as basis functions for approximate subspaces. In general, if one desires a higher order convergence for Galerkin-type approximations, choosing higher order basis functions gives higher convergence rate [33].,


    3. Application to nonlinear population balance equation

    In aerosol physics and environmental sciences, studying the flocculation of particles is widespread. The process of flocculation involves disperse particles in suspension combining into aggregates (i.e., a floc) and separating. The mathematical model used to study flocculation process is the well-known population balance equation (PBE) which describes the time-evolution of the particle size number density. The equations for the flocculation model track the time-evolution of the particle size number density u(t,x) and can be written as

    ut=F(u) (16)

    where

    F(u):=G(u)+A(u)+B(u),

    G denotes growth

    G(u):=x(gu)μ(x)u(t,x), (17)

    A denotes aggregation

    A(u):=12x0ka(xy,y)u(t,xy)u(t,y)dyu(t,x)ˉxx0ka(x,y)u(t,y)dy, (18)

    and B denotes breakage

    B(u):=¯xxΓ(x;y)kf(y)u(t,y)dy12kf(x)u(t,x). (19)

    The boundary condition is traditionally defined at the smallest size 0 and the initial condition is defined at t=0

    g(0)u(t,0)=¯x0q(x)u(t,x)dx,u(0,x)=u0(x)L1(0,¯x),

    where the renewal rate q(x) represents the number of new flocs entering the population. A floc is assumed to have a maximum size ¯x<. The function g(x) represents the average growth rate of the flocs of size x due to proliferation, and the coefficient μ(x) represents a size-dependent removal rate due to gravitational sedimentation and death. The function ka(x,y) is the aggregation kernel, which describes the rate with which the flocs of size x and y agglomerate to form a floc of size x+y. The fragmentation kernel kf(x) calculates the rate with which a floc of size x fragments. The integrable function Γ(x;y) represents the post-fragmentation probability density of daughter flocs for the fragmentation of the parent flocs of size y. In other words, all the fractions of daughter flocs formed upon the fragmentation of a parent floc sum to unity,

    y0Γ(x;y)dx=1 for all y(0,¯x]. (20)

    The population balance equation, presented in (16), is a generalization of many mathematical models appearing in the size-structured population modeling literature and has been widely used, e.g., to model the formation of clouds and smog in meteorology [49], the kinetics of polymerization in biochemistry [56], the clustering of planets, stars and galaxies in astrophysics [39], and even schooling of fish in marine sciences [47]. For example, when the fragmentation kernel is omitted, kf0, the flocculation model reduces to algal aggregation model used to describe evolution of a phytoplankton community [2]. When the removal and renewal rates are set to zero, the flocculation model simplifies to a model used to describe the proliferation of Klebsiella pneumoniae in a bloodstream [9]. Furthermore, the flocculation model, with only growth and fragmentation terms, was used to investigate the elongation of prion polymers in infected cells [17,26,18].

    The equation (16) has also been the focus of considerable mathematical analysis. Well-posedness of the general flocculation model was first established by Ackleh and Fitzpatrick [1,2] in an L2-space setting and later by Banasiak and Lamb [5] in an L1-space setting. Moreover, asymptotic behavior of the equation (16) has been a challenging task because of the nonlinearity introduced by the aggregation terms. Nevertheless, under suitable conditions on the kernels, the existence of a positive steady state has been established for the pure aggregation and fragmentation case [38]. For a review of further mathematical results, we refer readers to review articles by Menon and Pego [41], and Wattis [54] and the book by Ramkrishna [50]. Lastly, although the population balance equation has received substantial theoretical work, the derivation of analytical solutions for many realistic aggregation kernels has proven elusive. Towards this end, many discretization schemes for numerical simulations of the PBEs have been proposed. For instance, to approximate steady state solutions of PBEs, numerical schemes based on the least squares spectral method [23,24,25] and the finite element method [45,46,31] have been developed. For the further review of approximation methods we refer interested readers to the review by Bortz [8].


    3.1. Numerical implementation and results

    For the numerical implementation we adopt the scheme developed in Section 2.1. Therefore, the approximate formulation of (16) becomes the following system of n nonlinear ODEs for un=(α1,,αn)TRn :

    ˙un=Fn(un)=Gnun+PnA(Enun)+PnB(Enun), (21)
    un(0,x)=Pnu0(x), (22)

    where the matrix Gn is defined as in Section 2.3,

    PnA(Enun)=(α1n1j=1ka(xn1,xnj)αjΔx12ka(xn1,xn1)α1α1Δxα2n2j=1ka(xn2,xnj)αjΔx12n2j=1ka(xnj,xnn1j)αjαn1jΔxαn1ka(xnn1,xn1)α1Δx12n1j=1ka(xnj,xnnj)αjαnjΔx)

    and

    PnB(Enun)=(nj=2Γ(xn1;xnj)kf(xnj)αjΔx12kf(xn1)α1nj=3Γ(xn2;xnj)kf(xnj)αjΔx12kf(xn2)α2Γ(xnn1;xnn)kf(xnn)αnΔx12kf(xnn1)αn112kf(xnn)αn).

    The convergence of the approximate scheme (21)-(22) has been established in [1]. Therefore, the stationary solutions of the microbial flocculation model (16) can be systematically approximated by the stationary solutions of the system of nonlinear ODEs given in (21). We used Powell's hybrid root finding method [48] as implemented in Python 2.7.101 to find zeros of the steady state equation (see available code at [42]). For faster convergence rate, we provided the solver with the exact Jacobian of Fn(un) (see Section 2.2, Eqn 23) for the formulation of the Jacobian). For the purpose of illustration, for a post-fragmentation density function we chose the well-known Beta distribution 2 with α=β=2,

    Γ(x,y)=1[0,y](x)6x(yx)y3,

    1 scipy.optimize.fsolve

    2 Although normal and log-normal distributions are mostly used in the literature, Byrne et al. [16] have provided evidence that the Beta density function describes the fragmentation of small bacterial flocs.

    where 1I is the indicator function on the interval I. The aggregation kernel was chosen to describe flow within laminar shear field (i.e., orthokinetic aggregation)

    ka(x,y)=(x1/3+y1/3)3

    Other model rates were chosen arbitrarily as

    q(x)=a(x+1),g(x)=b(x+1)μ(x)=cxkf(x)=x,

    where a,b and c are some positive real numbers.

    The main advantage of this approximation scheme (21)-(22) is that it can be initialized very fast using Toeplitz matrices [40]. Fast initialization of the discretization scheme allows one to check the existence of the steady states at many discrete points efficiently. This in turn allows for the generation of the existence and stability regions of the steady states of the PBE in (16). To illustrate the existence regions of the steady states of the PBE, we discretized the intervals a[0,15],b[0,1] and c[0,5] with Δa=Δb=Δc=0.1. We note that for faster convergence the root finding method needs an initial seed close to the steady state solution. Since we have no information about the existing steady state, we seed the root finding method with 10 different uniform initial guesses i.e.,

    {u0(x)=2i|i=0,1,,9},

    before we conclude a positive steady state does not exist for a given point (a,b,c). Consequently, we checked for the existence of a positive steady state at each of these discrete points. As depicted in Figure 2a, approximate existence region of positive steady states of the PBE forms a three dimensional wedge like region. Moreover, in Figures 2b-2d, to deduce stability of each steady state solution, we checked the spectrum of the Jacobian matrix evaluated at each steady state. Particularly, if the real part of rightmost eigenvalue of the Jacobian matrix is negative, the steady state is identified as locally stable (blue region). Conversely, if the real part of the rightmost eigenvalue of the Jacobian matrix is positive the steady state is identified as unstable (red region). One can observe that growth (b) and removal (c) rates can balance the smaller renewal rates (a), and thus locally stable steady states exist. However, as the renewal rate gets larger steady states first become unstable and then cease to exist (yellow region). This is also illustrated in Figure 3b, where steady states start diverging for the larger renewal rates (a).

    Figure 2. Existence and stability regions for the steady states of the PBE a) Existence region for the steady states of the PBE forms a wedge like shape. b) Stability region for b=0.1, a[0,15] and c[0,5]. c) Stability region for b=0.5, a[0,15] and c[0,5]. d) Stability region for b=1.0, a[0,15] and c[0,5]. Color bar represents the real part of rightmost eigenvalue of the Jacobian matrix evaluated at each steady state. Yellow regions represents the region for which a positive steady state does not exists.
    Figure 3. An example steady-state solution of the PBE for b=0.5,a=c=1. b) Steady states for increasing renewal rate and b=c=1.

    Figure 3a illustrates an example stationary solution for b=0.5,a=c=1. To confirm that the function depicted in Figure 4 is indeed a locally stable steady state, we simulated the system of ODEs in (21)-(22) for t[0,10] with a collection of arbitrary initial conditions (Figure 4a) close to the steady state solution. One can observe in Figure 4 that the stationary solution is indeed locally stable and thus initial conditions, Figure 4, converge to the steady state depicted in Figure 4b. As depicted in Figure 4c and Figure 4d, convergence is also reflected in the evolution of the total number of flocs (zeroth moment),

    Figure 4. Time evolution of the flocculation model with arbitrary initial conditions. a) Four different initial conditions are chosen close to the steady state. b) Solution of the PBE for those initial conditions at t=10. c) Evolution of the total number M0(t) of the flocs for t[0,10]. d) Evolution of the total mass M1(t) of the flocs for t[0,10].
    M0(t)=ˉx0u(t,x)dxni=1xnixni1αiβni(x)dx=Δxni=1αi,

    and total mass of the flocs (first moment),

    M1(t)=ˉx0xu(t,x)dxni=1xnixni1αixβni(x)dx=Δx2ni=1αi(xni+xni1).

    Moreover, to confirm that the steady state solution is not changing with increasing dimension of approximate spaces Xn, we simulated our numerical scheme for different values of n. Figure 5a illustrates that stationary solutions converge to exact stationary solution of (16) as n. Furthermore, one can observe, in Figure (5) b that difference between approximate steady states for different values of n is considerably small.

    Figure 5. Change in zeroth and first moments with increasing dimension of the approximate space Xn. a) Change in the total number and the total mass of the flocs with respect to increasing dimension n. Dashed red lines and dotted green lines corresponds to the total number and the total mass of the flocs of the steady state for n=1000, respectively. b) Steady state solution for n=100 and n=500.

    3.2. Conditions for numerical stability of positive steady states

    In this section, we derive conditions for approximate local stability of the stationary solution of the nonlinear population balance equation defined in (16). In particular, we impose conditions on the model rates of the population balance equation for which the first statement of Corollary 2.2 holds. Towards this end, we use the well-known Gershgorin theorem for locating eigenvalues of a matrix. The Gershgorin theorem states that each eigenvalue of A lies in one of the the circular areas, called Gershgorin disks, in the complex plane that is centered at the ith diagonal element and whose radius is Ri,

    {zC:|zaii|Ri},

    where

    Ri=nj=1,ji|aji|.

    Since the approximate system for the microbial flocculation model is nonlinear, we linearize the system around its stationary solutions. Let uL1(0,¯x) be a stationary solution of (16) and denote the projection of the stationary solution u onto Rn by α=Pnu=[α1,,αn]T, then the Jacobian of the approximate operator Fn defined in (21) can be written as

    JF(α)=Gn+JA(α)+JB(α), (23)

    where Gn is defined in (13),

    JA(α)=(α1ka(xn1,xn1)Δxα1ka(xn1,xn2)Δxα1ka(xn1,xnn1)Δx0α2ka(xn2,xn1)Δxα2ka(xn2,xnn2)Δx00αN1ka(xnn1,xn1)Δx0000000)+(n1j=1ka(xn1,xnj)αjΔx000α1ka(xn1,xn1)Δxn2j=1ka(xn2,xnj)αjΔx00α2ka(xn1,xn2)Δxα1ka(xn2,xn1)Δx01j=1ka(xnn1,xnj)αjΔx0αn1ka(xn1,xnn1)Δxαn2ka(xn2,xnn2)Δxα1ka(xnn1,xn1)Δx0),

    and

    JB(α)=(12kf(xn1)Γ(xn1;xn2)kf(xn2)ΔxΓ(xn1;xn3)kf(xn3)ΔxΓ(xn1;xnn)kf(xnn)Δx012kf(xn2)Γ(xn2;xn3)kf(xn3)ΔxΓ(xn2;xnn)kf(xnn)Δx00012kf(xnn1)Γ(xnn1;xnn)kf(xnn)Δx00012kf(xnn))

    To bound the eigenvalues of JF(α) we use Gershgorin theorem. Consequently, the centers and the radii of Gershgorin disks are given by

    aii=1Δxg(xni)μ(xni)12kf(xni)αika(xni,xni)Δxnij=1ka(xni,xnj)αjΔx

    and

    Ri1Δxg(xni)+q(xni)+i1j=1Γ(xnj;xni)kf(xni)Δx+nij=1αjka(xnj,xni)Δx
    +nij=1,jiαjka(xni,xnj)Δx,

    respectively. Consequently, if we can show that

    |aii|>Rifor each i{1,,n}, (24)

    then each of the Gershgorin disks lie strictly on the left side of the complex plane. To this end, inequality (24) can be simplified as

    μ(xni)+12kf(xni)>q(xni)+i1j=1Γ(xnj;xni)kf(xni)Δx+nij=1αjka(xnj,xni)Δx (25)

    for each i{1,,n}. Accordingly, taking the limit of (25) as n yields

    μ(x)+12kf(x)>q(x)+x0Γ(y,x)kf(x)dy+¯xx0ka(x,y)u(y)dy (26)

    for all x[0,¯x] and together with the number conservation requirement (20) implies

    q(x)+12kf(x)μ(x)+¯xx0ka(x,y)u(y)dy<0

    for all x[0,¯x]. Conversely, note that the integral approximations in (25) are right Riemann sums. Therefore, if the functions Γ(y,x) and ka(x,y)u(y) are decreasing in y then integral approximations in (25) are under-approximations of the integrals in (26). Thus, the inequality stated in (26) ensures that the eigenvalues of the Jacobian JF(α) are strictly negative for all sufficiently large n. Now, we are in a position to summarize the results of this section in the following proposition.

    Proposition 3.1. Let u be a stationary solution of the microbial flocculation model (16). Moreover, if

    q(x)+12kf(x)μ(x)+¯xx0ka(x,y)u(y)dy<0 (27)

    for all x[0,¯x] and

    y(ka(x,y)u(y))0 and yΓ(y,x)0 (28)

    for all x[0,¯x] and y[0,¯x], then stationary solution u is approximately locally stable in the sense of Definition 2.1.

    To illustrate the utility of the framework developed in this section we applied our approach to the model rates given in Section 3.1 for generation of Figure 3a. The Beta distribution used for the post-fragmentation function Γ is not monotonically decreasing, and thus it does not satisfy the conditions of Proposition 3.1. However, Figure 6a-c illustrates that the model rates satisfy the conditions of Corollary 2.2. In fact, Figure 6d depicts that for the steady state illustrated in Figure 3a the eigenvalues of JF(Pnu) have strictly negative real part for n5. Therefore, as it has already been established in Figure 4, this steady state solution is locally asymptotically stable in the sense of Corollary 2.2.

    Figure 6. Eigenvalues of the Jacobian JF(α) multiplied by Δx for the steady state illustrated in Figure 3a. a) Eigenvalues of the Jacobian plotted in the complex plane for n=20. b) Eigenvalues of the Jacobian plotted in the complex plane for n=50. c) Eigenvalues of the Jacobian plotted in the complex plane for n=200. d) Change in the rightmost eigenvalue for increasing n. Dashed black line corresponds to the rightmost eigenvalue of the Jacobian for n=1000.

    4. Concluding remarks

    Our primary motivation in this paper was to show that available numerical tools in the literature can facilitate theoretical analysis of evolution equations. Towards this end we developed a numerical framework for theoretical analysis of evolution equations arising in population dynamical models. The main idea behind this framework is to approximate generators of semigroups of evolution equations and use numerical tools to study stability of steady states of evolution equations. We demonstrated the utility of our approach by applying the numerical framework to both linear and nonlinear size-structured population models. In particular, we generated approximate existence and stability regions of the steady states of both models (which can be difficult to obtain by using conventional analytical tools). We showed that our numerical framework can also be used to gain insight about the approximate local stability (see Definition 2.1) of stationary solutions. Furthermore, code used for the numerical simulations can be obtained from our GitHub repository under the open source MIT License (MIT) [42].

    Although the stability inequality in (7) holds for arbitrarily large finite time intervals, behavior of the solutions as t is unclear. Hence, we note that the local stability of the stationary solutions specified in Corollary 2.2 is not in a usual Lyapunov sense. In order to prove Lyapunov stability of stationary solutions using the approximation scheme presented in Section 2.1, one has to prove uniform convergence of the approximation scheme for all t0. Hence, as a future plan we wish to utilize the numerical framework presented here to establish Lyapunov stability of stationary solutions of general evolution equations.


    Acknowledgments

    Funding for this research was supported in part by grants NSF-DMS 1225878 and NIH-NIGMS 2R01GM069438-06A2. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS-0821794) and the University of Colorado Boulder. The Janus supercomputer is a joint effort of the University of Colorado Boulder, the University of Colorado Denver and the National Center for Atmospheric Research.


    [1] [ A.S. Ackleh, Parameter estimation in a structured algal coagulation-fragmentation model, Nonlinear Anal. Theory Methods Appl., 28 (1997): 837-854.
    [2] [ A.S. Ackleh,B.G. Fitzpatrick, Modeling aggregation and growth processes in an algal population model: analysis and computations, J. Math. Biol., 35 (1997): 480-502.
    [3] [ V. I. Arnold, Ordinary Differential Equations, Second printing of the 1992 edition. Universitext. Springer-Verlag, Berlin, 2006.
    [4] [ J. Banasiak, Blow-up of solutions to some coagulation and fragmentation equations with growth, Dyn. Syst., 1 (2011): 126-134.
    [5] [ J. Banasiak,W. Lamb, Coagulation, fragmentation and growth processes in a size structured population, Discrete Contin. Dyn. Syst. -Ser. B, 11 (2009): 563-585.
    [6] [ H.T. Banks,F. Kappel, Transformation semigroups andL 1-approximation for size structured population models, Semigroup Forum, 38 (1989): 141-155.
    [7] [ C. Biggs,P. Lant, Modelling activated sludge flocculation using population balances, Powder Technol., 124 (2002): 201-211.
    [8] [ D. M. Bortz, Chapter 17: Modeling and simulation for nanomaterials in fluids: Nanoparticle self-assembly. In Tewary, V. and Zhang, Y., editors, Modeling, characterization, and production of nanomaterials: Electronics, Photonics and Energy Applications, volume 73 of Woodhead Publishing Series in Electronic and Optical Materials, (2015), 419–441. Woodhead Publishing Ltd., Cambridge, UK.
    [9] [ D.M. Bortz,T.L. Jackson,K.A. Taylor,A.P. Thompson,J.G. Younger, Klebsiella pneumoniae Flocculation Dynamics, Bull. Math. Biol., 70 (2008): 745-768.
    [10] [ D. Breda, Solution operator approximations for characteristic roots of delay differential equations, Applied Numerical Mathematics, 56 (2006): 305-317.
    [11] [ D. Breda,O. Diekmann,S. Maset,R. Vermiglio, A numerical approach for investigating the stability of equilibria for structured population models, J. Biol. Dyn., 7 (2013): 4-20.
    [12] [ D. Breda,S. Maset,R. Vermiglio, Computing the characteristic roots for delay differential equations, IMA J Numer Anal, 24 (2004): 1-19.
    [13] [ D. Breda,S. Maset,R. Vermiglio, Pseudospectral Differencing Methods for Characteristic Roots of Delay Differential Equations, SIAM J. Sci. Comput., 27 (2005): 482-495.
    [14] [ D. Breda,S. Maset,R. Vermiglio, Pseudospectral approximation of eigenvalues of derivative operators with non-local boundary conditions, Applied Numerical Mathematics, 56 (2006): 318-331.
    [15] [ D. Breda,S. Maset,R. Vermiglio, Numerical approximation of characteristic values of partial retarded functional differential equations, Numer. Math., 113 (2009): 181-242.
    [16] [ E. Byrne, S. Dzul, M. Solomon, J. Younger and D. M. Bortz, Postfragmentation density function for bacterial aggregates in laminar flow, Phys. Rev. E, 83 (2011).
    [17] [ V. Calvez,M. Doumic,P. Gabriel, Self-similarity in a general aggregation-fragmentation problem, Application to fitness analysis, Journal de Math{é}matiques Pures et Appliqu{é}es, 98 (2012): 1-27.
    [18] [ V. Calvez,N. Lenuzza,M. Doumic,J.-P. Deslys,F. Mouthon,B. Perthame, Prion dynamics with size dependency-strain phenomena, J. Biol. Dyn., 4 (2010): 28-42.
    [19] [ A.M. De Roos, Demographic analysis of continuous-time life-history models, Ecol. Lett., 11 (2008): 1-15.
    [20] [ A. M. De Roos, PSPManalysis, 2014, https://staff.fnwi.uva.nl/a.m.deroos/PSPManalysis/index.html.
    [21] [ A.M. de Roos,O. Diekmann,P. Getto,M.A. Kirkilionis, Numerical equilibrium analysis for structured consumer resource models, Bulletin of Mathematical Biology, 72 (2010): 259-297.
    [22] [ O. Diekmann,M. Gyllenberg,J. A.J. Metz, Steady-state analysis of structured population models, Theoretical Population Biology, 63 (2003): 309-338.
    [23] [ C.A. Dorao,H.A. Jakobsen, Application of the least-squares method for solving population balance problems in Rd+1, Chemical Engineering Science, 61 (2006a): 5070-5081.
    [24] [ C.A. Dorao,H.A. Jakobsen, A least squares method for the solution of population balance problems, Computers & Chemical Engineering, 30 (2006b): 535-547.
    [25] [ C.A. Dorao,H.A. Jakobsen, Least-squares spectral method for solving advective population balance problems, Journal of Computational and Applied Mathematics, 201 (2007): 247-257.
    [26] [ M. Doumic-Jauffret and P. Gabriel, Eigenelements of a general aggregation-fragmentation model, arXiv: 0907.5467, 2009.
    [27] [ K. Engelborghs,T. Luzyanina,D. Roose, Numerical Bifurcation Analysis of Delay Differential Equations Using DDE-BIFTOOL, ACM Trans Math Softw, 28 (2002): 1-21.
    [28] [ J.Z. Farkas,T. Hagen, Stability and regularity results for a size-structured population model, J. Math. Anal. Appl., 328 (2007): 119-136.
    [29] [ J.Z. Farkas,P. Hinow, Steady states in hierarchical structured populations with distributed states at birth, Discrete Contin. Dyn. Syst. -Ser. B, 17 (2012): 2671-2689.
    [30] [ S. A. Gourley, R. Liu and J. Wu, Spatiotemporal Patterns of Disease Spread: Interaction of Physiological Structure, Spatial Movements, Disease Progression and Human Intervention, In Magal, P. and Ruan, S., editors, Structured Population Models in Biology and Epidemiology, number 1936 in Lecture Notes in Mathematics, (2008), 165–208. Springer Berlin Heidelberg.
    [31] [ M.J. Hounslow, A discretized population balance for continuous systems at steady state, AIChE J., 36 (1990): 106-116.
    [32] [ K. Ito,F. Kappel, The Trotter-Kato theorem and approximation of PDEs, Math Comp, 67 (1998): 21-44.
    [33] [ F. Kappel,K. Kunisch, Spline Approximations for Neutral Functional Differential Equations, SIAM J. Numer. Anal., 18 (1981): 1058-1080.
    [34] [ N. Kato, A Principle of Linearized Stability for Nonlinear Evolution Equations, Transactions of the American Mathematical Society, 347 (1995): 2851-2868.
    [35] [ T. Kato, Remarks on pseudo-resolvents and infinitesimal generators of semi-groups, Proc. Japan Acad., 35 (1959): 467-468.
    [36] [ T. Kato, Perturbation Theory for Linear Operators, Classics in Mathematics. Springer Berlin Heidelberg, Berlin, Heidelberg, 1976.
    [37] [ M.A. Kirkilionis,O. Diekmann,B. Lisser,M. Nool,B. Sommeijer,A.M. De Roos, Numerical continuation of equilibria of physiologically structured population models Ⅰ: Theory, Math. Models Methods Appl. Sci., 11 (2001): 1101-1127.
    [38] [ P. Laurencot,C. Walker, Steady states for a coagulation-fragmentation equation with volume scattering, SIAM Journal on Mathematical Analysis, 37 (2005): 531-548.
    [39] [ J. Makino,T. Fukushige,Y. Funato,E. Kokubo, On the mass distribution of planetesimals in the early runaway stage, New Astronomy, 3 (1998): 411-417.
    [40] [ S.A. Matveev,A.P. Smirnov,E.E. Tyrtyshnikov, A fast numerical method for the Cauchy problem for the Smoluchowski equation, Journal of Computational Physics, 282 (2015): 23-32.
    [41] [ G. Menon,R.L. Pego, Dynamical scaling in Smoluchowski's coagulation equations: Uniform convergence, SIAM Rev., 48 (2006): 745-768.
    [42] [ I. Mirzaev, Steady state approximation, 2015. https://github.com/MathBioCU/SteadyStateApproximation.
    [43] [ I. Mirzaev and D. M. Bortz, Criteria for linearized stability for a size-structured population model, arXiv: 1502.02754, 2015.
    [44] [ I. Mirzaev and D. M. Bortz, Stability of steady states for a class of flocculation equations with growth and removal, arXiv: 1507.07127, 2015 (submitted).
    [45] [ M. Nicmanis,M.J. Hounslow, Finite-element methods for steady-state population balance equations, AIChE J., 44 (1998): 2258-2272.
    [46] [ M. Nicmanis,M.J. Hounslow, Error estimation and control for the steady state population balance equation: 1. An a posteriori error estimate, Chemical Engineering Science, 57 (2002): 2253-2264.
    [47] [ H.-S. Niwa, School size statistics of fish, J. Theor. Biol., 195 (1998): 351-361.
    [48] [ M. Powell, A hybrid method for nonlinear equations, Numerical Methods for Nonlinear Algebraic Equations, null (1970): 87-114, {Gordon & Breach}.
    [49] [ H. R. Pruppacher and J. D. Klett, Microphysics of Clouds and Precipitation: Reprinted 1980, Springer Science & Business Media, 2012.
    [50] [ D. Ramkrishna, Population Balances: Theory and Applications to Particulate Systems in Engineering, Academic Press, 2000.
    [51] [ S.J. Schreiber,M.E. Ryan, Invasion speeds for structured populations in fluctuating environments, Theor. Ecol., 4 (2011): 423-434.
    [52] [ J.W. Sinko,W. Streifer, A New Model For Age-Size Structure of a Population, Ecology, 48 (1967): 910-918.
    [53] [ H.F. Trotter, Approximation of semi-groups of operators, Pacific J. Math., 8 (1958): 887-919.
    [54] [ J.A. Wattis, An introduction to mathematical models of coagulation-fragmentation processes: A discrete deterministic mean-field approach, Phys. Nonlinear Phenom., 222 (2006): 1-20.
    [55] [ G. F. Webb, Theory of Nonlinear Age-Dependent Population Dynamics, CRC Press, 1985.
    [56] [ R.M. Ziff,G. Stell, Kinetics of polymer gelation, J. Chem. Phys., 73 (1980): 3492-3499.
  • This article has been cited by:

    1. E. P. KIGHTLEY, A. PEARSON, J. A. EVANS, D. M. BORTZ, Fragmentation of biofilm-seeded bacterial aggregates in shear flow, 2018, 29, 0956-7925, 1062, 10.1017/S0956792518000074
  • Reader Comments
  • © 2017 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5202) PDF downloads(835) Cited by(1)

Article outline

Figures and Tables

Figures(6)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog