---
title: "PI Curves: Jassby and Platt"
author: "AJ Smit"
format:
  html:
    anchor-sections: true
    number-sections: true
    page-layout: article
---

::: callout-tip
## This Theory Accompanies the Following Lab

- [Lab 3. PI Curves](Lab3_PI_curves.qmd)
:::

# The Hyperbolic Tangent Model

The hyperbolic tangent model was proposed by Jassby and Platt (1976). It has become one of the most widely used models for describing the relationship between photosynthetic rate and irradiance (light intensity) in aquatic photosynthetic organisms, including algae ranging from kelp to phytoplankton. The model captures the core dynamics of photosynthesis, in which the rate of photosynthesis initially increases with light intensity but eventually saturates as the photosynthetic machinery reaches its maximum efficiency. This is a simple model, but it effective because the biologically meaningful parameters can be directly interpreted to assess plant or algal productivity in various light environments.

The hyperbolic tangent model is expressed as:

$$ P(I) = P_{\text{max}} \times \tanh\left(\frac{\alpha I}{P_{\text{max}}}\right) $$

Where:

- $P(I)$ represents the photosynthetic rate at a given irradiance $I$ (light intensity),
- $P_{\text{max}}$ is the maximum photosynthetic rate (also referred to as the light-saturated rate),
- $\alpha$ is the initial slope of the curve, which reflects the photosynthetic efficiency at low light levels,
- $I$ is the irradiance (light intensity),

The hyperbolic tangent function $\tanh$ is used to smoothly describe the transition between the linear increase in photosynthesis at low light intensities and the eventual plateau at higher intensities, where photosynthesis becomes light-saturated. The light compensation point, the point at which photosynthesis equals respiration (i.e., net photosynthesis is zero), can also be derived from this model.

The model describes the essential processes of photosynthesis with just two parameters: $P_{\text{max}}$ and $\alpha$. Both parameters are biologically meaningful and tell us how efficiently an organism can convert light into chemical energy under different light conditions. For example, higher values of $P_{\text{max}}$ indicate a greater potential for photosynthesis under optimal light conditions, while the value of $\alpha$ indicates how quickly photosynthesis responds to low light.

Applications of the hyperbolic tangent model are numerous. It is commonly used to estimate the photosynthetic performance of marine and freshwater algae, seagrasses, and macroalgae under varying environmental conditions. In kelp forests, for instance, we may use this model to assess how different species adapt to light intensities at various depths or how photosynthetic performance shifts in response to seasonal changes in light availability. Looking at phytoplankton, the model helps estimate productivity across different layers of the water column, where light intensity decreases with depth.

Below are a few lines of data taken from a hypothetical P-I experiment. The data are for five replicate experiments with the same light intensities (independent variable), representing conditions typically encountered by kelp at latitudes between -36° and -23°S.

```{r}
#| echo: false
#| messages: false
#| warning: false

# Load necessary libraries
library(minpack.lm) # For non-linear least squares fitting
library(ggplot2) # For plotting
library(knitr)

# Set seed for reproducibility
set.seed(666)

# Define parameters based on typical kelp conditions
Pmax_true <- 13 # True maximum photosynthetic rate (e.g., mg C m^-2 h^-1)
alpha_true <- 0.05 # True initial slope (photosynthetic efficiency)

# Light intensities typically encountered by kelp (μmol photons m^-2 s^-1)
light_intensities <- seq(0, 550, by = 50)

# Number of replicates
n_replicates <- 5

# Generate data for each replicate
data_list <- lapply(1:n_replicates, function(rep) {
  # Simulate true photosynthetic rates using the hyperbolic tangent model
  P_true <- Pmax_true * tanh(alpha_true * light_intensities / Pmax_true)
  
  # Add random noise to simulate experimental variability
  noise <- rnorm(length(P_true), mean=0, sd=0.5)  # Adjust sd for noise level
  P_observed <- round(P_true + noise, 2)
  
  # Create a data frame for the replicate
  data.frame(
    Replicate = rep,
    Light = light_intensities,
    Photosynthesis = P_observed
  )
})

# Combine all replicates into one data frame
data_all <- do.call(rbind, data_list)

# View 9 randomly selected rows from the data
random_rows <- data_all[sample(nrow(data_all), 9), ]

kable(random_rows, format = "markdown", digits = 2,
      col.names = c("Replicate", "Light (μmol photons m⁻² s⁻¹)", "Photosynthesis (mg C m⁻² h⁻¹)"))
```

```{r}
#| echo: false
#| messages: false
#| warning: false

# Fit the hyperbolic tangent model to each replicate
fit_results <- lapply(1:n_replicates, function(rep) {
  # Subset data for the replicate
  data_rep <- subset(data_all, Replicate == rep)
  
  # Define the model function
  model_formula <- Photosynthesis ~ Pmax * tanh(alpha * Light / Pmax)
  
  # Provide initial parameter estimates
  start_params <- list(Pmax = 10, alpha = 0.03)
  
  # Fit the model using non-linear least squares
  fit <- nlsLM(
    formula = model_formula,
    data = data_rep,
    start = start_params,
    control = nls.lm.control(maxiter = 1000)
  )
  
  # Return the fit object
  return(fit)
})
```

```{r}
#| echo: false
#| messages: false
#| warning: false
#| results: hide

# Summarise fit results for each replicate
for (i in 1:n_replicates) {
  cat("\n--- Replicate", i, "Fit Summary ---\n")
  print(summary(fit_results[[i]]))
}
```

```{r}
#| echo: false
#| messages: false
#| warning: false
#| label: fig-jassby_platt
#| fig.cap: "Nonlinear regression Jassby and Platt (1976) model fitted to simulated P-I data for a hypothetical kelp."

# Plot the data and fitted curves for all replicates
ggplot(data_all, aes(x = Light, y = Photosynthesis,
                     color = factor(Replicate))) +
  geom_point() +
  stat_function(
    fun = function(x) {
      coef(fit_results[[1]])['Pmax'] *
        tanh(coef(fit_results[[1]])['alpha'] *
               x / coef(fit_results[[1]])['Pmax'])
    },
    color = "black", size = 1, linetype = "dashed"
  ) +
  labs(
    title = "Photosynthesis-Irradiance Curves for Kelp (All Replicates)",
    x = "Light Intensity (μmol photons m⁻² s⁻¹)",
    y = "Photosynthetic Rate (mg C m⁻² h⁻¹)",
    color = "Replicate"
  ) +
  theme_minimal()
```

After fitting the model to the data, we can determine the values for $P_{\text{max}}$ and $\alpha$ for each replicate and determine the average value across the five fits. The combined plot (@fig-jassby_platt) displays the observed data points for all replicates and the fitted curve from the first replicate.

The average model fit values of the estimated parameters across all replicates are as follows:

- **$P_{\text{max}}$:** `r round(mean(sapply(fit_results, function(fit) coef(fit)['Pmax'])), 2)` mg C m⁻² h⁻¹
- **$\alpha$:** `r round(mean(sapply(fit_results, function(fit) coef(fit)['alpha'])), 2)` μmol photons m⁻² s⁻¹

# Considering the Light Compensation Point

The **light compensation point** ($I_c$) is the irradiance level at which the rate of photosynthesis equals the rate of respiration, resulting in a net photosynthetic rate of zero. Below this point, the organism consumes more energy (via respiration) than it produces through photosynthesis, leading to a net loss of energy. Estimating $I_c$ is important for determining the minimum light intensity required for the survival of photosynthetic organisms, after compensation for the effect of cellular respiration.

In the context of the **Jassby and Platt hyperbolic tangent model**, $I_c$ can be estimated by solving for the irradiance $I$ when the net photosynthetic rate $P(I)$ equals zero:

$$
0 = P_{\text{max}} \times \tanh\left(\frac{\alpha I_{\text{LCP}}}{P_{\text{max}}}\right)
$$

Since $\tanh(0) = 0$, the net photosynthetic rate is zero when $I = 0$. However, due to respiration, the net photosynthesis can be negative at zero light intensity. To account for respiration, we can modify the model to include **dark respiration rate** ($R$):

$$
P(I) = P_{\text{max}} \times \tanh\left(\frac{\alpha I}{P_{\text{max}}}\right) - R
$$

Now, $I_c$ is the irradiance at which $P(I) = 0$:

$$
0 = P_{\text{max}} \times \tanh\left(\frac{\alpha I_{\text{LCP}}}{P_{\text{max}}}\right) - R
$$

We can solve this equation numerically to find $I_{\text{LCP}}$.

```{r}
#| echo: false
#| messages: false
#| warning: false

# Set seed for reproducibility
set.seed(13)

# Define parameters based on typical kelp conditions
Pmax_true <- 13      # True maximum photosynthetic rate (e.g., mg C m^-2 h^-1)
alpha_true <- 0.05   # True initial slope (photosynthetic efficiency)
R_true <- 2          # True respiration rate (mg C m^-2 h^-1)

# Light intensities including finer increments around expected LCP
light_intensities <- seq(0, 550, by = 50)

# Number of replicates
n_replicates <- 5

# Generate data for each replicate
data_list <- lapply(1:n_replicates, function(rep) {
  # Simulate true net photosynthetic rates using the modified hyperbolic tangent model
  P_true <- Pmax_true * tanh(alpha_true * light_intensities / Pmax_true) - R_true
  
  # Add random noise to simulate experimental variability
  noise <- rnorm(length(P_true), mean = 0, sd = 0.5)  # Adjust sd for noise level
  P_observed <- P_true + noise
  
  # Create a data frame for the replicate
  data.frame(
    Replicate = rep,
    Light = light_intensities,
    Photosynthesis = P_observed
  )
})

# Combine all replicates into one data frame
data_all <- do.call(rbind, data_list)

# View 9 randomly selected rows from the data
random_rows <- data_all[sample(nrow(data_all), 9), ]

# Display the random rows in a nicely formatted table using kable
kable(random_rows, format = "markdown", digits = 2,
      col.names = c("Replicate", "Light (μmol photons m⁻² s⁻¹)", "Photosynthesis (mg C m⁻² h⁻¹)"))
```

```{r}
#| echo: false
#| messages: false
#| warning: false

# Fit the modified hyperbolic tangent model to each replicate and estimate LCP
fit_results2 <- lapply(1:n_replicates, function(rep) {
  # Subset data for the replicate
  data_rep <- subset(data_all, Replicate == rep)
  
  # Define the modified model function including respiration
  model_formula <- Photosynthesis ~ Pmax * tanh(alpha * Light / Pmax) - R
  
  # Provide initial parameter estimates
  start_params <- list(Pmax = 10, alpha = 0.03, R = 1)
  
  # Fit the model using non-linear least squares
  fit <- nlsLM(
    formula = model_formula,
    data = data_rep,
    start = start_params,
    lower = c(0, 0, 0),  # Parameters should be non-negative
    control = nls.lm.control(maxiter = 1000)
  )
  
  # Estimate the Light Compensation Point (LCP)
  # Solve for I when Photosynthesis = 0
  coef_fit <- coef(fit)
  Pmax_fit <- coef_fit['Pmax']
  alpha_fit <- coef_fit['alpha']
  R_fit <- coef_fit['R']
  
  # Define the function to find the root (Photosynthesis = 0)
  lcp_function <- function(I) Pmax_fit * tanh(alpha_fit * I / Pmax_fit) - R_fit
  
  # Use uniroot to find the LCP
  LCP <- uniroot(lcp_function, interval = c(0, 200))$root
  
  # Add LCP to the fit object
  fit$LCP <- LCP
  
  # Return the fit object
  return(fit)
})
```

```{r}
#| echo: false
#| messages: false
#| warning: false
#| results: hide

# Summarise fit results for each replicate
for (i in 1:n_replicates) {
  cat("\n--- Replicate", i, "Fit Summary ---\n")
  fit <- fit_results2[[i]]
  print(summary(fit))
  cat(sprintf("Estimated Light Compensation Point (LCP): %.2f μmol photons m⁻² s⁻¹\n", fit$LCP))
}
```

```{r}
#| echo: false
#| messages: false
#| warning: false
#| label: fig-jassby_platt2
#| fig.cap: "Nonlinear regression Jassby and Platt (1976) model fitted to simulated P-I data for a hypothetical kelp. This model includes the effect of the light compensation point."

# Plot the data, fitted curves, and LCP for all replicates
ggplot(data_all, aes(x = Light, y = Photosynthesis, color = factor(Replicate))) +
  geom_point() +
  stat_function(
    fun = function(x) {
      Pmax <- coef(fit_results2[[1]])['Pmax']
      alpha <- coef(fit_results2[[1]])['alpha']
      R <- coef(fit_results2[[1]])['R']
      Pmax * tanh(alpha * x / Pmax) - R
    },
    color = "black", size = 1, linetype = "dashed"
  ) +
  geom_vline(
    xintercept = fit_results2[[1]]$LCP,
    color = "red", linetype = "dotdash",
    show.legend = TRUE
  ) +
  geom_hline(yintercept = 0,
             color = "red", linetype = "dotted") +
  labs(
    title = "Photosynthesis-Irradiance Curves with Light Compensation Point",
    x = "Light Intensity (μmol photons m⁻² s⁻¹)",
    y = "Net Photosynthetic Rate (mg C m⁻² h⁻¹)",
    color = "Replicate"
  ) +
  annotate("text", x = fit_results2[[1]]$LCP, y = min(data_all$Photosynthesis) + 5,
           label = paste0("LCP = ", round(mean(sapply(fit_results2, function(fit) fit$LCP)), 2)),
           angle = 90, vjust = -0.5, color = "red") +
  theme_minimal()
```

The model fit to the data is in @fig-jassby_platt2. The average model fit values of the estimated parameters across all replicates are as follows:

- **$P_{\text{max}}$:** `r round(mean(sapply(fit_results2, function(fit) coef(fit)['Pmax'])), 2)` mg C m⁻² h⁻¹
- **$\alpha$:** `r round(mean(sapply(fit_results2, function(fit) coef(fit)['alpha'])), 2)` μmol photons m⁻² s⁻¹
- **$I_c$:** `r round(mean(sapply(fit_results2, function(fit) fit$LCP)), 2)` μmol photons m⁻² s⁻¹

# Platt et al. (1980) Model with Photoinhibition

Let's now look at the **Platt et al. (1980)** model, which incorporates photoinhibition into the photosynthesis-irradiance (P-I) relationship. This model extends the understanding of photosynthesis by accounting for the decrease in photosynthetic efficiency at high light intensities due to photoinhibition—a phenomenon where excessive light damages the photosynthetic apparatus, leading to reduced photosynthetic rates.

The model is expressed mathematically as:

$$
P(I) = P_{\text{max}} \left(1 - \exp\left(-\frac{\alpha I}{P_{\text{max}}}\right)\right) \exp\left(-\frac{\beta I}{P_{\text{max}}}\right)
$$

Where:

- $P_{\text{max}}$ is the maximum photosynthetic rate in the absence of photoinhibition.
- $\beta$ is the photoinhibition parameter (rate of decrease in photosynthesis at high light).
- $\exp$ denotes the exponential function.

This model combines the positive effect of light on photosynthesis at low irradiance with the negative effect of photoinhibition at high irradiance, providing a comprehensive description of the photosynthetic response across a wide range of light intensities.

```{r}
#| echo: false
#| messages: false
#| warning: false

set.seed(7)

# Define parameters based on typical kelp conditions
Pmax_true <- 13      # True maximum photosynthetic rate (e.g., mg C m^-2 h^-1)
alpha_true <- 0.05   # True initial slope (photosynthetic efficiency)
beta_true <- 0.0005  # True photoinhibition parameter
R_true <- 2          # True respiration rate (mg C m^-2 h^-1)

# Light intensities ranging to capture photoinhibition effects (μmol photons m^-2 s^-1)
light_intensities <- seq(0, 1600, by = 100)

# Number of replicates
n_replicates <- 5

# Generate data for each replicate
data_list <- lapply(1:n_replicates, function(rep) {
  # Simulate true net photosynthetic rates using the Platt et al. model
  P_true <- Pmax_true * (1 - exp(-alpha_true * light_intensities / Pmax_true)) *
            exp(-beta_true * light_intensities / Pmax_true) - R_true
  
  # Add random noise to simulate experimental variability
  noise <- rnorm(length(P_true), mean = 0, sd = 0.5)  # Adjust sd for noise level
  P_observed <- P_true + noise
  
  # Create a data frame for the replicate
  data.frame(
    Replicate = rep,
    Light = light_intensities,
    Photosynthesis = P_observed
  )
})

# Combine all replicates into one data frame
data_all <- do.call(rbind, data_list)

# View 9 randomly selected rows from the data
random_rows <- data_all[sample(nrow(data_all), 9), ]

# Display the random rows in a nicely formatted table using kable
kable(random_rows, format = "markdown", digits = 2,
      col.names = c("Replicate", "Light (μmol photons m⁻² s⁻¹)", "Photosynthesis (mg C m⁻² h⁻¹)"))
```


```{r}
#| echo: false
#| messages: false
#| warning: false

# Fit the Platt et al. model to each replicate
fit_results <- lapply(1:n_replicates, function(rep) {
  # Subset data for the replicate
  data_rep <- subset(data_all, Replicate == rep)
  
  # Define the model function including respiration
  model_formula <- Photosynthesis ~ Pmax * (1 - exp(-alpha * Light / Pmax)) *
                   exp(-beta * Light / Pmax) - R
  
  # Provide initial parameter estimates
  start_params <- list(Pmax = 10, alpha = 0.03, beta = 0.0001, R = 1)
  
  # Fit the model using non-linear least squares
  fit <- nlsLM(
    formula = model_formula,
    data = data_rep,
    start = start_params,
    lower = c(0, 0, 0, 0),  # Parameters should be non-negative
    control = nls.lm.control(maxiter = 1000)
  )
  
  # Return the fit object
  return(fit)
})
```

```{r}
#| echo: false
#| messages: false
#| warning: false
#| results: hide

# Summarise fit results for each replicate
for (i in 1:n_replicates) {
  cat("\n--- Replicate", i, "Fit Summary ---\n")
  fit <- fit_results[[i]]
  print(summary(fit))
}
```

```{r}
#| echo: false
#| messages: false
#| warning: false
#| label: fig-jassby
#| fig.cap: "Nonlinear regression Platt et al. (1980) model fitted to simulated P-I data for a hypothetical kelp. This model includes the effect of photoinhibition."

# Plot the data and fitted curves for all replicates
ggplot(data_all, aes(x = Light, y = Photosynthesis, color = factor(Replicate))) +
  geom_point() +
  stat_function(
    fun = function(x) {
      Pmax <- coef(fit_results[[1]])['Pmax']
      alpha <- coef(fit_results[[1]])['alpha']
      beta <- coef(fit_results[[1]])['beta']
      R <- coef(fit_results[[1]])['R']
      Pmax * (1 - exp(-alpha * x / Pmax)) * exp(-beta * x / Pmax) - R
    },
    color = "black", size = 1, linetype = "dashed"
  ) +
  geom_hline(yintercept = 0,
             color = "red", linetype = "dotted") +
  labs(
    title = "Photosynthesis-Irradiance Curves with Photoinhibition",
    x = "Light Intensity (μmol photons m⁻² s⁻¹)",
    y = "Net Photosynthetic Rate (mg C m⁻² h⁻¹)",
    color = "Replicate"
  ) +
  theme_minimal()
```




---

**References:**

- **Jassby, A. D., & Platt, T. (1976).** Mathematical formulation of the relationship between photosynthesis and light for phytoplankton. *Limnology and Oceanography*, 21(4), 540-547.
- **Platt, T., Gallegos, C. L., & Harrison, W. G. (1980).** Photoinhibition of photosynthesis in natural assemblages of marine phytoplankton. *Journal of Marine Research*, 38(4), 687-701.


