The Spatial Spead of Bark Beetle Infestations: Modelling and Simulation

Selwyn Hollis

Armstrong Atlantic State University

www.math.armstrong.edu/faculty/hollis

Abstract

We consider a reaction-diffusion-chemotaxis model for the spatial spread of bark-beetle infestations (specifically the southern pine beetle). Standard Keller-Siegel-type chemotaxis terms model the pheromonal attraction and repulsion of beetles to and from colonized trees. The hyperbolic character of these nonlinear chemotaxis terms results in very steep gradients and induces an amazing variety of complex behaviors, including patterns and traveling waves and pulses.

The mixed parabolic-hyperbolic character of the model makes numerical simulation a delicate problem. Following work of LeVeque et al., we apply a fractional-step method in which different methods specially suited for reaction, diffusion, and chemotaxis, respectively, are combined (without sacrificing second-order accuracy). We will discuss the implementation of this technique in Mathematica and present a number of generated animations.

Outline

Description of the model

Numerical Issues and methods

Mathematica implementation

Animations

Dendroctonus frontalis, the southern pine beetle (SPB).

Model

Variables

({{ φ, flying  beetles}, {ℓ, landed  beetles}, {𝒻, nesting beetles}, { A, aggregation pheromone}, { B, dispersal pheromone}, { ϑ, host tree odor}, {R, resin capacity}})

Assumptions

Forest Landscape

Attackable tree surface area: S(x, y)

Resin capacity: R(t, x, y) and (Subscript[R, 0])(x, y) = R(0, x, y)

R is an indicator of tree vigor, and tree death occurs as R->0 locally.

R is scaled such that (Subscript[R, max]) = 1 corresponds to complete resistance to attack.

Assumptions  (cont.)

Pheromones

Flying beetles ( φ ) disperse randomly and are chemotactically attracted by an aggregation pheromone ( A ), released by landed females ( ℓ ).

Flying beetles also are repelled by a “dispersal” pheromone ( B ) that is released by nesting males.

The sex ratio of nesting beetles is roughly 1:1; so there’s no need to track the male population.

We assume that A is produced at a rate proportional to ℓ, and B is produced at a rate proportional to 𝒻.

The net effect is a local change in φ at rate

∇·( φ ∇ψ(A, B) )

where we assume

\(TraditionalForm\`\(\(	\)\(ψ(A, 
        B) = β\/\(1 + B\) - \(\(α\/\(1 + 
                A\)\)\(.\)\)\)\)\)

Assumptions  (cont.)

Beetle-host interaction

In combination with the aggregation pheromone, tree odor ( ϑ )—i.e., volitiles released by tree resin—facilitates landing. As more holes are bored, more tree odor is released as the tree defends itself by filling holes with resin.

Landing rate of flying beetles: (Subscript[k, p]) S ϑ φ. 
(A female beetle is considered “landed” only if she is successful in drilling a hole.)

Landed beetles ( ℓ ) successfully attack (i.e., become nesting beetles) at rate (Subscript[k, f]) ϱ(R/(Subscript[R, 0])) ℓ and are killed by tree resin at rate (Subscript[k, f]) (1 - ϱ(
                  R/(Subscript[R, 0]))) ℓ .

Host odor volatiles ( ϑ ) are produced at rate ((Subscript[k, ϑ])1 
                    S + (Subscript[k, ϑ])2 ℓ) R, where 0 < (Subscript[k, ϑ])1 ≪ (Subscript[k, ϑ])2.

Resin capacity decreases due to attack at rate (Subscript[δ, R]) (ℓ + 𝒻) R and (slowly) recovers or declines at rate

(Subscript[k, R])(R/(Subscript[R, 0])) (1 - (R/(Subscript[R, 0])))((R/(Subscript[R, 0]))-ρ) ,           where 0 < ρ  < 1.

Assumptions  (cont.)

Emergence

Rate of emergence of new adults:

Q(𝒻)=\(\(4  σ\ η\^3\ \(S\^3\) 𝒻\)\/\(3 \( η\^4\) S\^4 + 𝒻\^4\)\)             (scramble competition)

This is a function of the “attack density” f/S, having a maximum value of σ S  when f/S = η. Thus η is a threshold parameter and σ is the maximum number of new adults per unit surface area.

Reëmergence

Nesting beetles ( 𝒻 ) reëmerge at rate r 𝒻.

The Differential Equations

({{flying , φ, (Subscript[φ, t])-𝒹 Δ φ, =, ∇·( φ ∇ψ(A, B) )+σ S (Subscript[Q, η])(𝒻/S)+r 𝒻-(Subscript[k, ℓ]) S ϑ φ-(Subscript[δ, φ]) φ}, {landed , ℓ, (Subscript[ℓ, t]), =, \(\(k\_ℓ\) 
                    S\ ϑ\ φ - \(k\_𝒻\) ℓ\)}, {nesting  , 𝒻, (Subscript[𝒻,  t]), =, (Subscript[k, 𝒻]) ϱ(
                      R/(Subscript[R, 0])) ℓ - ((Subscript[δ, 𝒻]) + 
                        r) 𝒻}, {aggregation pheromone, A, (Subscript[A, t])-ν Δ A, =, (Subscript[k, A]) ℓ - (Subscript[δ, A]) A}, {dispersal pheromone, B, (Subscript[B, t])-ν Δ B, =, (Subscript[k, B]) 𝒻 - (Subscript[δ, B]) B}, {host odor, ϑ, (Subscript[ϑ, t])-\(ν\ Δ ϑ\), =, \(\((k\_ϑ1 + \(k\_ϑ2\) ℓ)\)\ R - δ\_ϑ\ ϑ\)}, {resin capacity, R, (Subscript[R, t]), =, (Subscript[k, R])(R/(Subscript[R, 0])) (1-(R/(Subscript[R, 0])))((R/(Subscript[R, 0]))-ρ)-(Subscript[δ, R])(ℓ + 𝒻) R}})

Spatial domain:  [0, 1]×[0, 1]

Boundary condition for all diffusing species:  ∂(u/∂)n+γ u=0

Numerics

Issues

Diffusion equations are typically solved using by implicit methods because of severe stability constraints inherent in explicit methods.

The chemotaxis terms are nonlinear advection terms and have a hyperbolic character. So they are best treated using an (explicit) method for hyperbolic equations.

In the context of a typical Crank-Nicolson-type method applied to reaction-diffusion equations, steep gradients can cause spurious oscillations due to slow decay of high-frequency modes. With nonlinear chemotaxis, this problem can be severe.

Fractional-step Method  ( operator splitting )

View the system as

\(\(\(		\)\(\(∂\/∂t\) 
        Y = \((ℛ  + 𝒟  + 𝒜  )\) 
        Y\)\)\)

Naïve Splitting

Advance approximate solution from time (Subscript[t, n]) to time (Subscript[t, n + 1]) = (Subscript[t, n]) + h according to

\(\(\(		\)\(\(Y\_*\) = ℛ \((h)\) 
        Y\_n\)\)\)

\(\(\(		\)\(\(Y\_ ** \) = 𝒟 \((h)\) \(Y\_*\)\)\)\)

\(\(\(		\)\(Y\_\(n + 1\) = 𝒜 \((h)\) \(Y\_ ** \)\)\)\)

That is,

\(\(\(		\)\(Y\_\(n + 1\) = ℛ \((h)\) 𝒟 \((h)\) 𝒜 \((h)\) Y\_n\)\)\)

Strang Splitting

Second-order accuracy is preserved by a ”Strang splitting”:

\(\(\(		\)\(Y\_\(n + 1\) = ℛ \((h/
            2)\)\ \ 𝒟 \((h/
            2)\)\ \ 𝒜 \((h)\)\ \ 𝒟 \((h/
            2)\)\ \ ℛ \((h/2)\) Y\_n\)\)\)

Our present approach

Use TR-BDF2 for the diffusion step. (Trapezoid Rule + 2-step Backward Difference Formula)

For stepsize h,

{({{({{(Subscript[u, n + 1/2]), =, ((ℳ(ρ/
                                    4))^(-1))ℳ(-ρ/
                              4) (Subscript[u, n]) }, {(Subscript[u, n + 1]), =, ℳ(((ρ/3))^(-1)) ((4/3) (Subscript[u, n + 1/2])-(1/3) (Subscript[u, n]))}})}})

where

\(\(\(		\ \ \ \ \ \)\(ρ = \(\(h\/δ\^2\) 𝒹\ \ \ and\ \ ℳ \((μ)\) = ℐ - μ\ Δ\)\)\)\)

That is,

𝒟 (h)=ℳ(((ρ/3))^(-1)) ((4/3)ℳ(((ρ/4))^(-1))ℳ(-ρ/4) -(1/3) ℐ)

This method is “L-stable” as well as “A-stable,” allowing better control of high-frequency instability.

Use “Improved Euler” (TR) for the reaction step.

\(\(\(		\)\(Y\_\(n + 1\) = 
      Y\_n + \(h\/2\) \((ℛ \((Y\_n)\) + ℛ \((Y\_n + h\ ℛ \((Y\_n)\))\))\)\)\)\)

Use the MacCormack method for the chemotaxis step.

\(\(\(	\)\(\(φ\^*\) = φ - \(h\/δ\^2\) \((\((μ\_\(x\ f\)\ φ)\) \((∂\_\(x\ f\)ψ)\) - φ\ ∂\_\(x\ c\)ψ + \((μ\_\(y\ f\)\ φ)\) \((∂\_\(y\ f\)ψ)\) - φ\ ∂\_\(y\ c\)ψ)\)

    	\(φ\_+\) = \(1\/2\) \((\(φ\^*\) + φ)\) - \(h\/\(2  δ\^2\)\) \((\(φ\^*\)\ ∂\_\(x\ c\)ψ - \((μ\_\(x\ b\)\ \(φ\^*\))\) ψ + \(φ\^*\)\ ∂\_\(y\ c\)ψ - \((μ\_\(y\ b\)\ \(φ\^*\))\) \((∂\_\(y\ b\)ψ)\))\)\)\)\)

Diffusion step

ℳ[n_Integer, ρ_] := 
        SparseArray@
          Join[
Table[{i, i} -> 
                1 + (4 + γ) ρ, {i, 2, 
                n - 1}], 
Table[{i, i} -> 
                1 + (4 + γ) ρ, {i, (n^2) - n + 2, 
                (n^2) - 1}], 
Table[{i, i} -> 
                1 + (4 + γ) ρ, {i, n + 1, (n^2) - 2  n + 1, 
                n}], 
Table[{i, i} -> 
                1 + (4 + γ) ρ, {i, 2  n, (n - 1) n, 
                n}], 
Table[{i, i} -> 
                1 + (4 + 2  γ) ρ, {i, 1, n, 
                n - 1}], 
Table[{i, i} -> 
                1 + (4 + 2  γ) ρ, {i, (n^2) - n + 1, (n^2), 
                n - 1}], 
Table[{i, i} -> 
                1 + 4  ρ, {i, 
                (n^2)}], 
Table[{i, i + 1} -> 
                Which[Mod[i, n] == 0, 0, 
                  Mod[i, n] == 1, -2 ρ, 
                  True, -ρ], {i, 
                (n^2) - 1}], 
Table[{i, i - 1} -> 
                Which[Mod[i, n] == 0, -2 ρ, 
                  Mod[i, n] == 1, 0, True, -ρ], {i, 2, 
                (n^2)}], 
Table[{i, 
                  i + n} -> -2 ρ, {i, 
                n}], 
Table[{i, 
                  i + n} -> -ρ, {i, n + 1, 
                (n^2) - n}], 
Table[{i, 
                  i - n} -> -ρ, {i, n + 1, 
                (n^2) - n}], 
Table[{i, 
                  i - n} -> -2 ρ, {i, (n^2) - n + 1, 
                (n^2)}] ];

ℰ𝓊ℓℯ𝓇[
        u_, ρ_] := ℰ[Ϟ[ρ]] . 
        u;  Ϟ[r_Real] := 
      Round[(10^6) r];

\(\(𝒟𝒯𝓇𝓅𝓏[
          u_, ρ_] := 
        Lsolve["TR", Ϟ[ρ/
                2]]@ℰ𝓊ℓℯ𝓇[
            u, ρ/2];\)
\)
𝒟ℬℱ𝒟[
          u1_, u2_, ρ_] := 
        Lsolve["BFD", Ϟ[ρ/3]]@(4  u2 - u1 )/
          3 ;

\(\(𝒟[
          u_, ρ_] := 𝒟ℬℱ𝒟[
          u, 𝒟𝒯𝓇𝓅𝓏[
            u, ρ/2], ρ];\)\)

Chemotaxis step

\(\(χ[u_, v_, θ_, i_] := 
        Module[{𝓊, dvx, dvy}, 
          If[OddQ[i], 

𝓊 = 
              u - θ \((\((μfx . u)\) \((𝒹fx . v)\) - 
                      u \((dvx = 𝒹cx . v)\) + \((μfy . 
                            u)\) \((𝒹fy . v)\) - 
                      u \((dvy = 𝒹cy . 
                              v)\))\); 
 .5 \((𝓊 + 
                  u - θ \((𝓊*
                          dvx - \((μbx . 𝓊)\) \((𝒹bx . 
                              v)\) + 𝓊*
                          dvy - \((μby . 𝓊)\) \((𝒹by . 
                              v)\))\))\), 

𝓊 = 
              u - θ \((u \((dvx = 𝒹cx . v)\) - \((μbx . 
                            u)\) \((𝒹bx . v)\) + 
                      u \((dvy = 𝒹cy . v)\) - \((μby . 
                            u)\) \((𝒹by . 
                            v)\))\); 
 .5 \((𝓊 + 
                  u - θ \((\((μfx . 𝓊)\) \((𝒹fx . v)\) - 𝓊*dvx + \((μfy . 𝓊)\) \((𝒹fy . 
                              v)\) - 𝓊*
                          dvy)\))\)\ \ ]\ \ ];\)\)

mμfx[n_Integer] := 
        SparseArray@
          Join[
Table[{i, i} ->  .5, {i, 1, 
                (n^2)}], 
Table[{i, i + 1} -> 
                Which[Mod[i, n] == 0, 0, True,  .5], {i, 
                (n^2) - 1}], 
Table[{i, 
                  i - 1} ->  .5, {i, n, (n^2), n}] ];
mμfy[n_Integer] := 
        SparseArray@
          Join[
Table[{i, i} ->  .5, {i, 1, 
                (n^2)}], 
Table[{i, i + n} ->  .5, {i, 
                1, (n^2) - n}], 
Table[{i, 
                  i - n} ->  .5, {i, (n^2) - n + 1, 
                (n^2)}] ];
mμbx[n_Integer] := 
        SparseArray@
          Join[
Table[{i, i} ->  .5, {i, 1, 
                (n^2)}], 
Table[{i, i - 1} -> 
                Which[Mod[i, n] == 1, 0, True,  .5], {i, 2, 
                (n^2)}], 
Table[{i, i + 1} ->  .5, {i, 
                1, (n^2) - n + 1, n}] ];
mμby[n_Integer] := 
        SparseArray@
          Join[
Table[{i, i} ->  .5, {i, 1, 
                (n^2)}], 
Table[{i, i - n} ->  .5, {i, 
                n + 1, (n^2)}], 
Table[{i, 
                  i + n} ->  .5, {i, 1, n}] ];
m𝒹fx[n_Integer] := 
        SparseArray@
          Join[
Table[{i, i} -> -1, {i, 1, 
                (n^2)}], 
Table[{i, i + 1} -> 
                Which[Mod[i, n] == 0, 0, True, 1], {i, 
                (n^2) - 1}], 
Table[{i, i - 1} -> 
                1, {i, n, (n^2), n}] ];
m𝒹fy[n_Integer] := 
        SparseArray@
          Join[
Table[{i, i} -> -1, {i, 1, 
                (n^2)}], 
Table[{i, i + n} -> 1, {i, 1, 
                (n^2) - n}], 
Table[{i, i - n} -> 
                1, {i, (n^2) - n + 1, (n^2)}] ];
m𝒹cx[n_Integer] := 
        SparseArray[
          Join[
Table[{i, i + 1} -> 
                Which[MemberQ[{0, 1}, Mod[i, n]], 0, True,  .5], {i, 2, 
                (n^2) - 1}], 
Table[{i, i - 1} -> 
                Which[MemberQ[{0, 1}, Mod[i, n]], 0, True, - .5], {i, 2, 
                (n^2)}] ], {(n^2), (n^2)}];
m𝒹cy[n_Integer] := 
        SparseArray[
          Join[
Table[{i, i + n} ->  .5, {i, n + 1, 
                (n^2) - n}], 
Table[{i, 
                  i - n} -> - .5, {i, n + 1, (n^2) - n}] ], {(n^2), 
            (n^2)}];
m𝒹bx[n_Integer] := 
        SparseArray@
          Join[
Table[{i, i} -> 1, {i, 1, 
                (n^2)}], 
Table[{i, i - 1} -> 
                Which[Mod[i, n] == 1, 0, True, -1], {i, 2, 
                (n^2)}], 
Table[{i, 
                  i + 1} -> -1, {i, 1, (n^2) - n + 1, 
                n}] ];
m𝒹by[n_Integer] := 
        SparseArray@
          Join[
Table[{i, i} -> 1, {i, 1, 
                (n^2)}], 
Table[{i, 
                  i - n} -> -1, {i, n + 1, 
                (n^2)}], 
Table[{i, 
                  i + n} -> -1, {i, 1, n}] ];

Setup of matrices and LinearSolve objects

\(setup[n_, 𝒹𝓉_, 𝒹_List] := 
      With[{ρ = 𝒹𝓉*\(n\^2\) 𝒹}, 

Map[\((Lsolve[
                  "TR", Ϟ@#] = 
                LinearSolve[ℳ[n, #]])\) &, ρ/
            8]; 

Map[\((Lsolve[
                  "BFD", Ϟ@#] = 
                LinearSolve[ℳ[n, #]])\) &, ρ/
            6]; 

Map[\((ℰ[Ϟ@\(Abs@#\)] = ℳ[n, #])\) &, \(-ρ\)/
            8]\ ; 

μfx = 
          mμfx[n]; μfy = mμfy[n]; μbx = 
          mμbx[n]; μby = 
          mμby[n]; 

𝒹fx = 
          mμfx[n]; 𝒹fy = m𝒹fy[n]; 𝒹bx = 
          m𝒹bx[
            n]; 

𝒹by = 
          mμby[n]; 𝒹cx = m𝒹cx[n]; 𝒹cy = 
          m𝒹cy[n];]\)

Clear[Lsolve,ℰ];
γ=2;
setup[n=100, 𝒹𝓉=.001, {𝒹=.002, (𝓂=50)𝒹}];

Setting parameters and compiling functions

α =  .45; β =  .45;
σ = 15; η =  .5;
δφ =  .5; kℓ = 20/5; δℓ = 20; 
    kf = δℓ; δf = 4;
ka = 200; δa = 20; kb = 200; δb = 20; 
    kϑ1 = 1; 
    kϑ2 = 200*5; δϑ = 20;
kr = 10; δr = 40;
\(\(Q = 
        Compile[{{𝒻, _Real, 
              1}}, \(4\ η\^3\ s\^3\ 𝒻\)\/\(3 \( η\^4\) s\^4 + 𝒻\^4\), {{s, _Real, 1}}];\)\)
ϱ = 
        Compile[{{x, _Real, 1}}, 
          (2. /π) ArcTan[((((x/(1 + (10. ^(-12)) - x))))^8)]];
Clear[φ, 𝒻, ϑ, ℓ, a, b, 
        R];
Vars = {φ, 𝒻, ϑ, ℓ, a, b, 
          R};
vars = {#, _Real, 1} & /@ Vars;
varnames = ToString /@ Vars;

    
\(\(\(F["χ"] = 
        Compile[{{a, _Real, 1}, {b, _Real, 
              1}}, \(-α\)\ \((1 + a)\)\^\(-1\) + β \((1 + b)\)\^\(-1\)];\)\(
\)
    \)\)
\(\(F["φ"] = 
        Compile[Evaluate@vars, 
          s\ σ\ Q[𝒻] + δf\ 𝒻/3 - 
            kℓ\ s\ ϑ\ φ - δφ\ φ, {{s, _Real, 1}}];\)\)
\(\(F["ℓ"] = 
        Compile[Evaluate@vars, 
          kℓ\ s\ ϑ\ φ - δℓ\ ℓ, {{s, _Real, 1}}];\)\)
F["\[ScriptF]"] = 
        Compile[Evaluate@vars, 
          kf ϱ[
                1 - R/R0] ℓ - δf 𝒻, {{R0, _Real, 
              1}}];
F["R"] = 
        Compile[Evaluate@vars, 
          kr  R (1 - 
                  R/R0) (R/
                    R0 -  .5) - δr (ℓ  + 𝒻) 
              R, {{R0, _Real, 1}} ];
\(\(F["ϑ"] = 
        Compile[Evaluate@vars, 
          kϑ1\ R + 
            kϑ2ℓ\ R\  - δϑ\ ϑ\ ];\)\)
F["a"] = 
        Compile[Evaluate@vars, ka ℓ - δa a];
F["b"] = 
        Compile[Evaluate@vars, 
          kb 𝒻 - δb b ];

Main loop

Do[
Do[
  
\(𝓇ℯ𝒶𝒸𝓉1 = \(Apply[Apply, #] &\) /@ \((\({F[#], Vars} &\) /@ varnames)\)\);
\(ℯ𝓊ℓℯ𝓇𝓈𝓉ℯ𝓅 = 
                    Vars +  .5  𝒹𝓉\ 𝓇ℯ𝒶𝒸𝓉1\);
\(𝓇ℯ𝒶𝒸𝓉2 = \(Apply[
                          Apply, #] &\) /@ \((\({F[#], ℯ𝓊ℓℯ𝓇𝓈𝓉ℯ𝓅} &\) /@ 
                          varnames)\)\);
\({φ, 𝒻, ϑ, ℓ, a, b, R} = 
                    Vars +  .25  𝒹𝓉 \((𝓇ℯ𝒶𝒸𝓉1 + 𝓇ℯ𝒶𝒸𝓉2)\)\);

  
\({ϑ, a, 
                      b} = \(𝒟[#,  .5  𝓂\ 𝒹*𝒹𝓉*n\^2] &\) /@ {ϑ, a, b}\);
\(φ = 𝒟[φ,  .5  𝒹*𝒹𝓉*n\^2]\);

  
\(φ = χ[φ, ψ = \(F["χ"]\)[a, b], 𝒹𝓉*n\^2, j]\);

  
\({ϑ, a, 
                      b} = \(𝒟[#,  .5  𝓂\ 𝒹*𝒹𝓉*n\^2] &\) /@ {ϑ, a, b}\);
\(φ = 𝒟[φ,  .5  𝒹*𝒹𝓉*n\^2]\);

  
\(𝓇ℯ𝒶𝒸𝓉1 = \(Apply[Apply, #] &\) /@ \((\({F[#], Vars} &\) /@ 
                          varnames)\)\);
\(ℯ𝓊ℓℯ𝓇𝓈𝓉ℯ𝓅 = 
                    Vars +  .5  𝒹𝓉\ 𝓇ℯ𝒶𝒸𝓉1\);
\(𝓇ℯ𝒶𝒸𝓉2 = \(Apply[
                          Apply, #] &\) /@ \((\({F[#], ℯ𝓊ℓℯ𝓇𝓈𝓉ℯ𝓅} &\) /@ 
                          varnames)\)\);
\({φ, 𝒻, ϑ, ℓ, a, b, R} = 
                    Vars +  .25  𝒹𝓉 \((𝓇ℯ𝒶𝒸𝓉1 + 𝓇ℯ𝒶𝒸𝓉2)\)\),

{j, 25}];

makeplot[],

{100}];

Selected References

Bark beetles

www.barkbeetles.org

Logan, White, Bentz, and Powell: Model Analysis of Spatial Patterns in Mountain Pine Beetle Outbreaks. Theoretical Population Biology 53 (1998)

Turchin, Lorio, Taylor, and Billings: Why do populations of southern pine beetles fluctuate? Environmental Ecology 20 (1991) See also: USDA Forest Service.

Reeve, Rhodes, and Turchin: Scramble competiton in the southern pine beetle, Dendroctonus frontalis. Environmental Ecology 23 (1998) See also: USDA Forest Service.

J. Hughes: Modeling the effect of landscape pattern on mountain pine beetles. Master’s thesis, Simon Fraser University (2002)

Numerics

Tyson, Stern, and LeVeque: Fractional-step methods applied to a chemotaxis model. J. Math. Biology 41 (2000)

R. LeVeque: Numerical Methods for Conservation Laws. Birkhäuser (1992)

G. Strang: On the construction and comparison of difference schemes. SIAM JNA 5 (1968)