https://rpubs.com/ScottWeiner19/1042972

Raw Data Plot

Here is a look at the raw data. It is 30 optical density readings an hour apart.

Source: https://www.jove.com/v/10511/growth-curves-generating-growth-curves-using-colony-forming-units

Bacterial growth is an example of logistic growth and takes place in 4 phases:

Lag, Exponential, Stationary, and Death.

library(readxl)
  #loads 'readxl' package
  
ODforR <- read.csv('ODforR_CSV.csv')
  #importing csv into R dataframe and naming datafram "ODforR" (read.csv is a function in the readxl library)

Hrs <-ODforR$`Time..hr.`
OD <- ODforR$OD600
plot(Hrs,OD, 
     main = expression(italic("E.coli \n")* plain("culture OD600 Growth Curve \n (with growth phases)")),
     xlab= "Time (hours)",
     ylab= "Optical Density (OD600)",
     pch=20)
     
     
#Next lines superimpose text onto plot


text(8, 0.2, "Lag", col="blue", cex=0.8)
text(15, 2.2, "Exponential", col="darkgreen", cex=0.8)
text(3,3.8, "Stationary", col="#FFA500", cex=0.8)
text(29, 2.9, "Death", col="red", cex=0.8)
  #first two numbers or "arguments" are xy coordinates for text

#next lines add arrows 

arrows(8, 0, 4, 0, length=0.15, col="blue")
arrows(14, 2, 9, 2, length=0.15, col="darkgreen")
arrows(3, 4, 7, 4, length=0.15, col="#FFA500")
arrows(29, 3.1, 29, 3.6, length=0.15, col="red")

  #first two numbers or "arguments" are xy coordinates for arrow start
  #second two arugments are xy coordiantes for arrow end
  #legnth sets length of the arrow head

Exponential For Loop Simulation

Linearizing raw data to estimate k and half life

#only looking at first 10 data points (before plateau)
logdatafit <- lm(log(OD[1:10])~Hrs[1:10])


#slope of the line will give estimate of K
k <- coef(logdatafit)[2]
t2 <- 0.693/k

#cat() function concatenates string(plain text) +  variable

cat("k = ",k)
## k =  0.7571047
cat("Doubling Time = ",t2)
## Doubling Time =  0.9153291

Initializing Rest of Parameters

#initializing time to 0 and starting OD
t <- 0
A <- 0.005


N = 60
Time = 10



dt <- Time/N
  #setting the time step size (interval between each     #discrete time point)
  #every iteration of loop represents 10 minutes



for(i in 2:N) {    
  t[i] <-  t[i-1] + dt
  A[i] <- A[i-1] + k*A[i-1]*dt     
}
  #for every iteration of this equation, where i will equal 2 going to N (60), values of t and A get updated
  #both t and A become variables storing a vector (of all of the values calculated for t and A)

print(t)
##  [1] 0.0000000 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333 1.0000000
##  [8] 1.1666667 1.3333333 1.5000000 1.6666667 1.8333333 2.0000000 2.1666667
## [15] 2.3333333 2.5000000 2.6666667 2.8333333 3.0000000 3.1666667 3.3333333
## [22] 3.5000000 3.6666667 3.8333333 4.0000000 4.1666667 4.3333333 4.5000000
## [29] 4.6666667 4.8333333 5.0000000 5.1666667 5.3333333 5.5000000 5.6666667
## [36] 5.8333333 6.0000000 6.1666667 6.3333333 6.5000000 6.6666667 6.8333333
## [43] 7.0000000 7.1666667 7.3333333 7.5000000 7.6666667 7.8333333 8.0000000
## [50] 8.1666667 8.3333333 8.5000000 8.6666667 8.8333333 9.0000000 9.1666667
## [57] 9.3333333 9.5000000 9.6666667 9.8333333
print(A)
##  [1] 0.005000000 0.005630921 0.006341453 0.007141644 0.008042806 0.009057681
##  [7] 0.010200616 0.011487772 0.012937346 0.014569834 0.016408316 0.018478785
## [13] 0.020810514 0.023436471 0.026393781 0.029724257 0.033474987 0.037698998
## [19] 0.042456014 0.047813288 0.053846566 0.060641148 0.068293098 0.076910603
## [25] 0.086615500 0.097545001 0.109853631 0.123715415 0.139326337 0.156907108
## [31] 0.176706294 0.199003823 0.224114946 0.252394694 0.284242897 0.320109838
## [37] 0.360502617 0.405992324 0.457222110 0.514916281 0.579890540 0.653063520
## [43] 0.735469767 0.828274375 0.932789451 1.050492670 1.183048167 1.332330063
## [49] 1.500448964 1.689781801 1.903005437 2.143134510 2.413564059 2.718117524
## [55] 3.061100801 3.447363122 3.882365616 4.372258518 4.923968126 5.545294728
plot(Hrs,OD,
     pch = 20,
     main= "Exponential Model Simulation",
     xlab= "Time (hours)",
     ylab= "Optical Density (OD600)",
     xlim = c(0,Time), 
     ylim = c(0,4), 
)

lines (t,A,
       col = "red")



Ac  <-  A[1] * exp(k*t)
  #Ac is the theoretical exponential solution using only starting OD,
  #estimated K, and the array of t values
  #calculated in the for loop above
print(Ac)
##  [1] 0.005000000 0.005672455 0.006435350 0.007300846 0.008282745 0.009396700
##  [7] 0.010660472 0.012094209 0.013720772 0.015566093 0.017659593 0.020034650
## [13] 0.022729131 0.025785995 0.029253980 0.033188378 0.037651918 0.042715763
## [19] 0.048460650 0.054978173 0.062372245 0.070760752 0.080277439 0.091074035
## [25] 0.103322677 0.117218650 0.132983508 0.150868598 0.171159072 0.194178432
## [31] 0.220293691 0.249921217 0.283533381 0.321666079 0.364927283 0.414006731
## [37] 0.469686924 0.532855605 0.604519907 0.685822416 0.778059382 0.882701394
## [43] 1.001416819 1.136098403 1.288893454 1.462238069 1.658895981 1.882002619
## [49] 2.135115100 2.422268941 2.748042398 3.117629466 3.536922681 4.012607074
## [55] 4.552266754 5.164505823 5.859085559 6.647080042 7.541052719 8.555256708
lines(t,Ac,
      col="blue",
      lwd= 2)

legend("topleft", legend = c("Raw Data", "Exponential Simulation", "Theoretical Solution"), 
       col = c("black","red","blue"), 
       lty = c(0,1,1), #line type
       pch=c(20,-1, -1), #point type
       cex = 1.2 #text size
)

Logistical Model For Loop

Now will use all 30 time points, incorporating carrying capacity as well as updating k estimate

Cap <-  4

t <- 0    

A <- 0.005    

N = 100

Time = 30  

dt <- Time/N  


k <-  1.2
  #this k was adjusted to get closest fit, visually


for(i in 2:N) 
  
{    t[i] <-  t[i-1] + dt


A[i] <- A[i-1] + k*A[i-1]* (1-A[i-1]/Cap) * dt       

}

  #raw data plot
plot(Hrs,OD,
     
     xlim = c(0,Time),
     ylim = c(0,4), 
     xlab = "Time(hours)",
     ylab= "Optical Density (OD600)",
     main ="Logistic Model Simulation"
)
  #logistic numerical model simulation 
lines(t,A,
      col="red",
      lty=2
)


Ac  <-  A[1] * exp(k*t)
  #theoretical solution

lines(t,Ac, col="blue")

legend("bottomright", legend = c("Raw Data", "Logistic Simulation", "Theoretical Solution"),
       col = c("black","red","blue"), 
       lty = c(0,2,1), 
       pch=c(20,-1, -1), 
       cex = 1.2
)

Non-Linear Least Squares “nls()”

Now we are using the nls function, to predict values using non-linear regression. The function runs an algorithm that attempts various values for the variables of the equation, and chooses the set of values that yields smallest difference between the predicted and observed data.

y <- OD

x <- Hrs
  #redefining x,y for easier reuse in future






tryfit <- nls(y ~ a / (1 + exp(-b * (x - c))), start = list(a = 1, b = 1, c = 1))

  # a = max value 
  # b = k value 
  # c = midpoint 

  #lists starting estimates for a,b,c


summary(tryfit)
## 
## Formula: y ~ a/(1 + exp(-b * (x - c)))
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a  3.87283    0.03201  120.98  < 2e-16 ***
## b  0.92742    0.07340   12.63 4.37e-13 ***
## c  6.55497    0.09774   67.07  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1431 on 28 degrees of freedom
## 
## Number of iterations to convergence: 12 
## Achieved convergence tolerance: 4.512e-06
pred <- predict(tryfit)

#plotting raw data

plot(Hrs,OD,
     xlim = c(0,Time),
     ylim = c(0,4), 
     xlab = "Time(hours)",
     ylab= "Optical Density (OD600)",
     main ="NLS Model vs Logistic Numerical Model",
     pch=20)


#Adding NLS model
lines(Hrs,pred,
      col="darkgreen",
      lwd=2,
      )

#Adding logistic numerical method model (Euler)
lines(t,A,
      col = "red",
      lty= 2
      )

legend("bottomright", legend = c("Raw Data", "NLS Model", "Logistic \nNumerical Model"), 
       col = c("black","darkgreen","red"), lty = c(0,1,2), pch=c(20,-1, -1), cex = 1.2)