Objectives

  1. To understand basic concepts of differential equation models.
  2. To generate & plot simulated individual time series using an ODE solver.
  3. To generate individual differences in simulated time series.

Outline

Introduction
1. Motivation for dynamic systems in developmental research
2. Brief introduction to differential equations
Modeling with Diff Eq 1. Constant growth (linear model) 2. Proportional growth/decay (exponential model) 3. Growth to an asymptote (logistic model)

Introduction: Dynamical systems theories and describing change

The study of change is an integral part of developmental science. Development or change occurs over long spans of time (decades), short time spans (seconds or less), and time scales in between. As developmental scientists, we theorize that developmental changes occur in more or less a lawful form, initiated, moderated, and/or regulated by forces within and outside of an individual.

In developmental theory, we hear a lot about “dynamical systems”. A “dynamic” system is characterized by constant change or progress. There can be “dynamics” in individual development, between two individuals, or between a group of people (e.g. a family system). Individual and group development can be measured on short or long time scales, depending on the phenomenon of interest. However, it may be overly vague to simply state that a system has “dynamics”. We need to be specific about how the system changes and how these interrelationships are defined. This specificity can be achieved by assigning a specific mathematical form to the nature of the change. This is where differential equations are applied to dynamical systems theory. In brief, a differential equation is a function that describes how a variable changes over a period of time relative to itself and/or other parameters. This is in contrast to traditional growth modeling, where the growth function describes the overall shape (or functional form) of the growth curve.

Introduction: What is a differential equation?

A differential equation is an expression which relates a function \(x(t)\) to its derivative(s). A derivative describes the instantaneous slope (the tangent) along any point of \(x(t)\). Derivatives are written with notation \[\frac{dx}{dt}=\frac{"Change \space in \space x"}{"Change \space in \space t"}\] As you can see, a differential equation looks a lot like the equation for a slope (\(m=\frac{x2-x1}{t2-t1}\). dx and dt are defined as infintesimal changes in x and t.

First-order differential equations (meaning that only the first derivative is involved) can be generally written as \[\frac{dx}{dt}=G(x(t))\] Where \(G(x(t))\) is some function of \(x(t)\).

For example, a simple differential equation is: \[\frac{dx}{dt}=10x(t)\] This is stating that the derivative of \(x(t)\) is equal to 10 times the original function \(x(t)\). \(x(t)\) is changing in proportion to itself at any given time.

We can find a solution to this differential equation by finding a function \(x(t)\) that satisfies the conditions of the differential equation (\(x(t)\) is the functional form of the growth). The process of solving differential equations is somewhat akin to solving algebraic problems (e.g. solving \(4+3~x=12\) for x), except that the solution is a function, not a number.

For example, one possible solution (among many) to \(\frac{dx}{dt}=10x(t)\) is \(x(t)=e^{10t}\). We can verify this is a solution by differentiating \(x(t)\) with respect to t: \[\frac{dx}{dt}=10~e^{10t}\] Which is \(10*x(t)\). Later, I will show the process of solving for the general solution of \(\frac{dx}{dt}=10x(t)\).

It is also useful to understand the form of growth in terms of the differential equation itself. We can interpret the differential equation as an expression of how a system is changing in time. That is, the differential equation is an explicit statement of the “rules” of change a system is following. In some cases, the differential form may be more useful to developmentalists than the functional form of growth: 1. When the differential form allows us to translate theories of change processes in systems into mathematical forms that are more obvious than functional forms. 2. When the differential form is simpler than the functional form (reduced order polynomials) 3. In some cases, the differential form is all that exits (i.e. there is no exact solution). In these cases, we can still approximate time series from the differential form.

In the following examples, we will only deal with ordinary differential equations (ODEs), or differential equations with only one independent variable. There are also partial differential equations (PDEs), which contain more than one possible independent variables, and stochastic differential equations (SDEs), which add an element of random noise to a dynamic process.

Modeling with diff eq examples:

In this section, we will describe three (simple) ordinary differential equation models. We will discuss how the models map onto possible theories of change processes.

We will also show how to simulate time series from the diff eq models using an ODE solver from the deSolve package. An ODE solver calculates the trajectory of a system over time, given an initial value, a differential equation, and a time span. The ODE solver calculates this trajectory through the process of numerical integration, an algorithm for approximating the next value of the system, given at least one prior value and the differential equation.

Example 1: Constant growth (Linear model)

This form of growth is the most simple, and gives us a great way to practice translating theory into math.

Suppose we have a theory that some variable changed in a constant way. For example, we could model the amount of money in a savings account, if $10 was added every day. The change in the account balance is always $10/day, regardless of past or current balances.

Thsi can be expressed in mathematical form:
\[\frac{dx}{dt}=10\] where we have defined dt to be one day.

More generally, constant growth is expressed as: \[\frac{dx}{dt}=a\] where a is some constant number. If a is positive, the person learns new words every day, if a is negative, the person “loses” words every day, and if a is zero, there is no change in vocabulary over time.

The properties of this model makes it easily recognizable as equivalent to linear growth. You can prove this equality by solving the differential equation with the technique of separable differential equations:
Multiply both sides by dt: \[dx=a~dt\] And integrate: \[\int_{-\infty}^{\infty} dx=\int_{-\infty}^{\infty} a~dt\] \[x(t)=a~t+x_0\] Where \(x_0\) is the intercept, or initial value of the system at \(t=0\).

Simluating time series from the constant change diff eq

Step one: Load libraries and define the differential equation in a function

#load libraries
library(deSolve)
## 
## Attaching package: 'deSolve'
## The following object is masked from 'package:graphics':
## 
##     matplot
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.3.1
#Define the exponential growth diff eq function
lineargrowth.ode <- function(time, x, parms){
  with(as.list(c(x,parms)),{
    dxdt <- a
    list(dxdt)
  })
}

Note that when we define the function in R, we include time in the arguments: this is required for the solver, but is not used in the actual equation.

Step two: Define the parameter value (growth rate), initial value, and times over which you wish to obtain a solution.
In this case we’ve defined the initial amount as $1, the growth parameter as 10 dollars/day, and we ask for the value of the system over 50 days.

#Define an initial value
initvalues <- c(x = 1)
#Define the model parameter: intrinsic growth rate
parms <- c(a = 10)
#Define the times over which you want a solution
times <- seq(0, 50, 1)

Step three: Plug these values and function into the lsoda() function to obtain a numerical solution

#use ode() to obtain solution
out <- lsoda(y = initvalues, times = times, func = lineargrowth.ode, parms = parms) 
out=as.data.frame(out)

Step four: Plot the output

ggplot(out, aes(x=time, y=x))+geom_point()+
  xlab("Time (days)")+
  ylab("money")+
  ggtitle("Example 1: Constant rate")

This is a fairly easy result to obtain algebraically, but the output from the ODE solver shows us the solution we would expect. Using an ODE solver for this problem may feel like overkill, but later examples will deal with less obvious solutions, and for certain advanced problems, the ODE solver may be the only way to simulate a time series.

Example 2: Proportional growth/decay (Exponential model)

This model is the simplest model that proposes that change in a system is relative to the variable of interest. The model of proportional growth/decay is written as: \[\frac{dx}{dt}=a~x\] Where a is some constant. If a is positive, the system grows proportionately to itself. If a is negative, the system decays proportionately to itself.

For example, an intuitive application of proportional growth comes from population studies. Suppose we begin with a population of rabbits. Assuming unlimited resources, the rabbits reproduce with some average litter size. This leads to new generations of rabbits, who also reproduce, meaning that as the population increases, the rate of growth also increases.

This might sound familiar as an exponential model of growth/decay. We can prove the equivalency by solving the differential equation: \[\frac{1}{x}~dx=a~dt\] Then integrate: \[\int_{} \frac{1}{x}dx=\int_{} a~dt\] \[ln(x)=at+C\] Where C is an arbitrary constant.
Then exponentiate both sides to get x by itself: \[e^{ln(x)}=e^{at+C}\] \[x=e^{at}e^C\] \[x(t)=x_0 e^{at}, \space (x_0=e^C)\]

Simulating time series from proportional change diff eq

Step 1:Define the exponential growth diff eq as a function

#Define the exponential growth diff eq function
expgrowth.ode <- function(time, x, parms){
  with(as.list(c(x,parms)),{
    dxdt <- a*x
    list(dxdt)
  })
}

Step 2: Define the growth parameter, initial value, and times over you wish to obtain a solution In this case, we assume the initial condition is 1 at time \(t=0\), and the rate of growth is defined by: \[\frac{dx}{dt}=0.2~x\]

#Define an initial value
initvalues <- c(x = 1)
#Define the model parameter: intrinsic growth rate
parms <- c(a = .2)
#Define the times over which you want a solution
times <- seq(0, 50, 1)

Step 3: Plug these values and function into the lsoda() function to obtain a numerical solution

#use ode() to obtain solution
out <- lsoda(y = initvalues, times = times, func = expgrowth.ode, parms = parms) 
out1=as.data.frame(out)

Step four: Plot the output

ggplot(out1, aes(x=time, y=x))+geom_point()+
  xlab("Time")+
  ylab("x")+
  ggtitle("Proportional growth")

Simulating trajectories with between person differences

Up to this point, we have simulated only individual time series. However, our ODE solver also allows us to simulate multiple time series with different parameters/starting values. This way, we can simulate between-person differences.

Example 1: Different growth parameters

#make a list of multiple values for the parameter
outlist <- list()
plist <- cbind(a = runif(10,min=.05,max=.2))
#check out what parameters were generated
plist
##                a
##  [1,] 0.12731073
##  [2,] 0.08751583
##  [3,] 0.16605328
##  [4,] 0.09486285
##  [5,] 0.13491786
##  [6,] 0.14215061
##  [7,] 0.14133901
##  [8,] 0.17523753
##  [9,] 0.08659222
## [10,] 0.11182122
#define initial value and times
initvalues.1=c(x=10)
times=seq(0,50,1)
#Loop through the matrix of parameter values
for(i in 1:nrow(plist)){
  outlist[[i]] <- lsoda(y = initvalues.1, times = times, func = expgrowth.ode, parms = plist[i,])
}
#plot multiple outputs
plot(out,outlist, lwd=2,main="trajectories with different growth parameters")

Example 2: Different initial values

#define multiple initial values
outlist.1=list()
val.list=cbind(x=runif(10,min=1,max=250))
#check out what initial values were generated
val.list
##                x
##  [1,] 193.561432
##  [2,] 113.538465
##  [3,]   3.926156
##  [4,] 157.824274
##  [5,] 136.439653
##  [6,] 143.294892
##  [7,]  54.639607
##  [8,] 213.666358
##  [9,] 180.471389
## [10,] 233.847350
#Now define times and parameter values
times=seq(0,10,.25)
parms=c(a=.2)
#run model
for(i in 1:nrow(val.list)){
  outlist.1[[i]] <- ode(y=val.list[i,],times=times,func=expgrowth.ode,parms=parms)
}
#plot output
plot(out,outlist.1, xlim=c(0,10), ylim=c(0,500), lwd=2,main="trajectories with different initial values")

Example 3: Growth to an asymptote (Logistic model)

For some processes, we may theororize that growth occurs in an accelerating manner at first (similar to exponential growth), but eventually progress slows until some “cap” is reached. For example, in population studies, population growth may occur proportionately at first, but as resources become limited, growth slows (e.g. smaller litter sizes, less reproduction/survival of offspring to maturity).

This model is expressed mathematically as: \[\frac{dx}{dt}=a~x~(1-\frac{x}{K})\] Where a is a positive growth parameter as before, and K is the upper limit (called the carrying capacity).

You can see how this equation is similar to exponential growth in early stages. If x is small relative to K, then \(\frac{x}{K}\) will be close to zero, meaning that \((1-\frac{x}{K})\) will be close to 1. The whole expression is therefore close to the equation for exponential growth. As x increases relative to K, the proportional growth of the system slows. Finally, when \(x=K\), the growth of the system is zero.

Simulating a single trajectory with logistic growth

Step 1: Define logistic ODE

#define the logistic function
logistic <- function(time, x, parms){
  with(as.list(c(x,parms)),{
    dx <- a * x * (1 - x / K)
    list(dx)
  })
}

Step 2: Define parameters, initial value, and times

#Initial Value
initvalues.1 <- c(x = 0.1)
#Model Parameters: intrinsic growth rate r and carrying capacity K
parms.1 <- c(a = 0.1, K = 10) 
#Times over which you want a solution
times <- seq(0, 100, 1)

Step 3: Run model

#Run model
out.2 <- ode(y = initvalues.1, times = times, func = logistic, parms = parms.1) 
out.2a=as.data.frame(out.2)

Step 4: Plot output

#plot output
ggplot(out.2a, aes(x=time, y=x))+geom_point()+
  xlab("Time")+
  ylab("x")+
  ggtitle("Logistic growth")

Plot multiple trajectories with logistic model

Different growth parameters, same carrying capacity

#generate multiple growth parameters, r, with carrying capacity K=10
outlist.2 <- list()
plist.2 <- cbind(a = runif(10,min=0.05,max=2),
               K = 10)
#check out what r's were generated
plist.2
##               a  K
##  [1,] 1.0439389 10
##  [2,] 1.9734983 10
##  [3,] 1.6113277 10
##  [4,] 1.9641446 10
##  [5,] 0.4608000 10
##  [6,] 0.9048871 10
##  [7,] 0.7292276 10
##  [8,] 1.6292653 10
##  [9,] 0.7128629 10
## [10,] 0.6938174 10
#Loop through the matrix of parameter values
for(i in 1:nrow(plist.2)){
  outlist.2[[i]] <- ode(y = initvalues.1, times = times, func = logistic, parms = plist.2[i,])
}
#plot output
plot(out.2,outlist.2, lwd=2)

Conclusions:

We have discussed several basic ODE models, and how they might be applicable to theories of change.

ODE solvers are useful tools for generating simulated data from ODE models. We have shown here that the ODE form produces the same time series as the functional form of several models.

In order to more realisitically simulate data, one would add a certain amount of random noise to the “ideal” data generated by the ODE solver. This will be demonstrated in future tutorials.