With the release of version 6.3 of the COMSOL Multiphysics® software, it is now possible to use previous solutions from a time-dependent study at the current time step. That is, the solution from an earlier time can be used to affect the solution at the current time. This has many applications, including for the modeling of delay differential equations. Let’s take a look at this new functionality in the context of a somewhat whimsical example.
Mountain Meadow Musings
A while ago, I had the chance to do a bit of hiking in the mountains and spotted a marmot. After taking a few snapshots, I started thinking about the life cycle of this adorable herbivore. My thoughts quickly drifted to coming up with an equation to describe the population change of marmots munching in the mountain meadows.
The fuzzy marmot that served as inspiration for this blog post.
I thought back to some of my college textbooks on numerical methods and remembered all manner of equations that could be useful, including the Lotka–Volterra equations and the logistic differential equation. Taking inspiration from these, let’s postulate that we should write two ordinary differential equations (ODEs): one for the rate of change of the number of marmots in the meadow, M(t), and another for the rate of change of the amount of edible biomass, B(t).
The first equation, for the rate of change in the marmot population, has terms related to reproduction, R(t), and to death, D(t), which we will look into in more detail shortly. The second equation has three constants: G_0 is the rate of growth, B_{max} is the maximum mass of vegetation that can grow in the meadow, and E_0 is the marmots’ rate of consumption per day.
Now, turning to the equation for the marmot population, we will start with a reproduction term considering a linear dependence upon the current population:
The marmot death rate, on the other hand, should be related to available resources, so as an intermediate step, we will define the forage area per day:
where A_0 is the meadow size. We will assume that the greater the forage area, the greater the chances of encountering a predator. We will also assume a fixed number of predators per area, P(t) = P_0. So, the marmot death rate is proportional to the current number of marmots, M(t), each marmot’s daily forage area, and the predator density. This leads to an equation for the death rate that depends upon both unknowns:
We will further assume that the entire colony of marmots comes out of hibernation at once, and we will solve these equations to predict the population change over the summertime, after which the colony nests back into its burrow for the winter. We will also include a stop condition to stop the solver if the number of marmots or the amount of vegetation drops below a critical value, implying colony or environmental collapse.
Screenshot showing how to model a set of ODEs. Different ODEs can have different units. A Stop Condition feature is also added to the Time-Dependent Solver feature.
These types of equations can be solved in COMSOL Multiphysics® using the Global ODEs and DAEs interface, as shown in the screenshot above. This interface defines the two variables that we are solving for, M(t) and B(t), as well as the equation for each in terms of their time derivatives and their initial values, M_\text{init}=M(t=0) and B_\text{init}=B(t=0). We can solve these equations to predict how the population evolves over the summer.
Solutions to the equations for the dynamics of the vegetation biomass (green) and marmot population (brown), with an annotation of the population at the end of the summer. The marmot population initially rises steeply and then falls due to increased predation.
The results show that the marmot population initially rises sharply and rapidly consumes the available food. This leads to a correction, a fall in the population due to increased predation. After some more time the population of marmots reaches an equilibrium. (As a side note for the interested reader: Determining if this is a stable equilibrium can be quite challenging.)
Delay Differential Equations Account for Additional Factors
At this point we have some plausible-looking results, but we might also have some concerns. There are at least two significant factors that should be considered:
- Marmots start mating after coming out of hibernation but do not reproduce instantaneously; there is a gestational period.
- If the marmot population increases for a sufficiently long time, then more predators will be attracted into the meadow.
Let’s look at how we might describe these two factors mathematically and how to implement them within our model. To address the first, we will want to make the reproduction rate a function of the population at the time of conception and consider the gestation period, T_g, during which the reproduction rate is zero. So, our equation for reproduction rate becomes:
where H(t-T_g) is the Heaviside step function. The rate is also reduced by the fraction \frac{M\left( T_g \right)}{M_\text{init}}, i.e., the total population decline over the gestation period. We’ll assume that once the young marmots start foraging, the predation of the older marmots ceases. Therefore, this fraction stays fixed after t=T_g.
We now have a delay differential equation that can be entered into the software using a combination of an if() statement with the at(time,variable) operator. We we can write the reproduction term as:
R0*if(t > Tg,at(t-Tg,M)*at(Tg,M)/Minit,0)
The first argument to the at() operator is the previous time at which we want to evaluate the second argument. To enable this functionality of referring to previous times during the solution, we need to modify the solver settings as shown in the screenshot below. Using Strict time stepping with an output time step no larger than the delay time ensures that we capture the start of the birthing season. Otherwise, the model is solved in the same way as before.
Screenshot showing how to enable the use of time operators during the solution.
Let’s take a look at the population when we incorporate this nonconstant reproduction rate. We observe first the population decline over the period when no marmots are born and then a rise in population that starts to level off.
Population of marmots (brown) and amount of vegetation (green) when a delay term is introduced into the birth rate.
Finally, let’s modify the death rate to account for the increase in predators that occurs if the average number of marmots increases for a sufficiently long time span. That is, we need to track the time-averaged marmot population over the previous few days.
Rather than evaluate this integral over all previous time steps over the period T_p, we can apply the fundamental theorem of calculus and write an equation for the rate of change of the average population over the previous days:
This equation could be solved simply by adding another global equation to our model, but it would be evaluated starting from t=0, so we will use another if() statement to define M(t<0) = M_\text{init}. We use the moving average to scale up the predator density P(t) = \left( P_0 /M_\text{init}\right) M_\text{ave}(t), thus affecting the death rate. Solving again, we see the effect on the population at the end of the summer.
Good news! The population has increased, and the marmots can burrow down for their long winter hibernation.
Results showing the marmot population (brown) and biomass amount (green) with both the reproduction rate and the death rate being dependent upon previous solutions.
Experimental observation confirms a rise in marmot population.
A Brief Remark About Population Balance Modeling
It should be emphasized at this point that the example we have put together here is meant to demonstrate some exciting new functionality in the context of an easy to understand model. In fact, real population dynamics models can be quite a bit more sophisticated than this, often using a compartmentalized or interval-based approach, where different variables are used to track different fractions of the total population. An introductory example to this approach is the SEIR model used in epidemiology (Susceptible, Exposed, Infectious, Recovered and immune), as discussed in our blog post “Modeling the Spread of COVID-19 with COMSOL Multiphysics®”.
This interval-based approach is particularly of interest in chemical engineering, where it is often used to track the sizes of droplets and precipitates. The following tutorial models highlight this type of modeling:
- Oil–Water Flow Through an Orifice — A Droplet Population Model
- Precipitation of Barium Sulfate
- Crystallization of Benzoic Acid in a Mixed Suspension, Mixed Product Removal Crystallizer
Caveats and Closing Remarks
We have shown the fundamentals of implementing and solving delay differential equations in COMSOL Multiphysics® version 6.3. Although this simple model can produce some very interesting dynamics, keep in mind that this example is meant for learning purposes. Once you understand this simple example, you will be ready to tackle more sophisticated modeling tasks. You might, for example, want to combine these delay differential equations in time with differential equations in space. Every interface in the software supports these new operators. Regardless of the combination of interfaces that you want to use in your multiphysics modeling, the technique for referring to previous solutions will be exactly the same as what we’ve shown here.
Want to try out the model discussed above? Download the related MPH-file in the Application Gallery:
Comments (0)