Hello everyone!

In this post I will show how to estimate fungicide \(EC_{50}\) values in R. In dose-response terms, \(EC_{50}\) is the concentration at which a compound inhibits 50% of growth relative to the untreated control.

In plant pathology, \(EC_{50}\) is widely used to evaluate fungal sensitivity to fungicides. An increase in \(EC_{50}\) for part of a population may indicate selection toward less sensitive isolates.

Load the packages

The core model-fitting engine used here is drc, which provides several dose-response models. For stratified, multi-isolate workflows, I also recommend the ec50estimator project, which was designed specifically to streamline \(EC_{50}\) estimation in plant pathology datasets. If you do not have the packages installed yet, run install.packages("drc") in your console and install the others as needed.

library(tidyverse)
library(drc)
library(ec50estimator)
library(cowplot)
library(ggridges)
library(ggthemes)

Data

First, we need a dataset from which \(EC_{50}\) can be estimated. Here I will simulate data using a standard fungicide dose-response relationship (Noel et al. 2018). If you already have your own dataset, you can import it directly in R and move to the modeling steps.

\[ f(x) = c +\frac{d-c}{1+exp(b(log(x)-log(e)))}\]

This is the four-parameter log-logistic model, LL.4(). Its parameters are b, c, d, and e. Parameter b is the slope around the inflection point, c is the lower asymptote, d is the upper asymptote, and e corresponds to the \(EC_{50}\). Variable x represents the fungicide dose or concentration.

Using this model, we can create a helper function to simulate mycelial growth values.

grow = function(dose, b,c,d,e){
  x = dose
  y = c + (d-c)/(1+exp(b*(log(x)-log(e))))
  erro = (y * 0.02 * ( (d + 0.05 * d) - y ))
 
  
  growth = rnorm(length(x),y, erro)
  return(growth)
}

Below is a simple example. I create a vector of doses called doses, simulate growth values, and inspect the resulting data frame.

# create dose vector
doses = rep(c(0, 1e-6,1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1),5)
set.seed(1) # keeps the simulated values reproducible

# simulate mycelial growth
growth = grow(doses, b = 1, c = 0, d = 20, e = 0.001)

df1 = data.frame(doses,growth)
head(df1)

  doses    growth
1 0e+00 19.749418
2 1e-06 20.054870
3 1e-05 19.405505
4 1e-04 19.816651
5 1e-03 10.724917
6 1e-02  1.245888

We can then visualize the simulated observations.

df1 %>% 
  ggplot()+
  geom_jitter(aes(factor(doses),growth, color = growth),
              width = 0.05,
              size = 3.5, 
              alpha = 1)+
  geom_boxplot(aes(factor(doses),growth), fill = NA, size = .7)+
  scale_color_viridis_c(option = "D",direction = 1)+
  labs(x ="Dose", y = "Mycelial growth")+
  theme_minimal_hgrid(font_size = 12)+
  theme(legend.position = "none")

Estimating \(EC_{50}\) for a single isolate

To estimate \(EC_{50}\), we first fit an appropriate dose-response model. Here we use drm() from the drc package. We provide:

the formula, which maps response and dose columns
the data frame
the model function used for fitting

In this example, the formula is growth ~ doses, the data frame is df1, and the model is LL.4().

model1 = drm(formula = growth~doses,  data = df1, fct = LL.4())
summary(model1)


Model fitted: Log-logistic (ED50 as parameter) (4 parms)

Parameter estimates:

                 Estimate  Std. Error  t-value   p-value    
b:(Intercept)  1.0447e+00  7.9231e-02  13.1857 2.415e-15 ***
c:(Intercept) -1.7187e-02  2.3108e-01  -0.0744    0.9411    
d:(Intercept)  2.0030e+01  1.8943e-01 105.7401 < 2.2e-16 ***
e:(Intercept)  9.7556e-04  6.5352e-05  14.9277 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.7130597 (36 degrees of freedom)

The fitted curve can be inspected visually.

plot(model1)

To obtain the \(EC_{50}\) estimate itself, we use the function ED(). The argument respLev = 50 requests the effective dose corresponding to 50% response reduction. Confidence intervals can also be obtained with interval = "delta". This low-level workflow is useful because it makes the underlying estimation step explicit before moving to larger automated pipelines such as ec50estimator.

ED(model1, respLev=50, interval = "delta")


Estimated effective doses

         Estimate Std. Error      Lower      Upper
e:1:50 9.7556e-04 6.5352e-05 8.4302e-04 1.1081e-03

Estimating \(EC_{50}\) for multiple isolates

In practice, we often need to estimate \(EC_{50}\) for many isolates rather than a single one. To illustrate this, I will create a hypothetical dataset with 50 fungal isolates, half collected in an organic field and half in a conventional field. These isolates will be evaluated against two fungicides, here called Fungicide A and Fungicide B.

Fungicide A

For Fungicide A, I assume isolates from the conventional field are less sensitive.

set.seed(1)
n_isolates = 25

nrep = 5
n_dose =7
n_field = 2
n = n_isolates*nrep*n_dose*n_field
field = c(rep("Organic",nrep*n_dose), rep("Conventional",nrep*n_dose))
isolate = sort(rep(1:((n_isolates*2)),nrep*n_dose))
doses = (rep(c(0,1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1),nrep*n_field))
EC50 = c(5e-3, 1e-1)


fungicide_A = data.frame(isolate,
           field = rep(field,n_isolates),
           fungicide = "Fungicide A",
           dose = rep(doses,n_isolates)) %>% 
  mutate(effect = case_when(
    field == "Conventional" ~ rnorm(n,EC50[2],0.0001),
    field == "Organic" ~ rnorm(n,EC50[1],0.0001))) %>% 
  mutate(growth = grow(dose, b = 1, c = 0, d = 20, e = effect))%>% 
  dplyr::select(-effect)

head(fungicide_A)

  isolate   field   fungicide  dose     growth
1       1 Organic Fungicide A 0e+00 20.2082399
2       1 Organic Fungicide A 1e-05 20.1168279
3       1 Organic Fungicide A 1e-04 19.2479678
4       1 Organic Fungicide A 1e-03 15.8123455
5       1 Organic Fungicide A 1e-02  7.3206757
6       1 Organic Fungicide A 1e-01  0.6985264

Fungicide B

For Fungicide B, I assume both field types have the same sensitivity pattern.

set.seed(1)
n_isolates = 25
nrep = 5
n_dose =7
n_field = 2
n = n_isolates*nrep*n_dose*n_field
field = c(rep("Organic",nrep*n_dose), rep("Conventional",nrep*n_dose))
isolate = sort(rep(1:((n_isolates*2)),nrep*n_dose))
doses = (rep(c(0,1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1),nrep*n_field))
EC50 = c(5e-3, 5e-3)


fungicide_B = data.frame(isolate,
                        field = rep(field,n_isolates),
                        fungicide = "Fungicide B",
                        dose = rep(doses,n_isolates)) %>% 
  mutate(effect = case_when(
    field == "Conventional" ~ rnorm(n,EC50[2],0.0001),
    field == "Organic" ~ rnorm(n,EC50[1],0.0001))) %>% 
  mutate(growth = grow(dose, b = 1, c = 0, d = 20, e = effect)) %>% 
  dplyr::select(-effect)

head(fungicide_B)

  isolate   field   fungicide  dose     growth
1       1 Organic Fungicide B 0e+00 20.2082399
2       1 Organic Fungicide B 1e-05 20.1168279
3       1 Organic Fungicide B 1e-04 19.2479678
4       1 Organic Fungicide B 1e-03 15.8123455
5       1 Organic Fungicide B 1e-02  7.3206757
6       1 Organic Fungicide B 1e-01  0.6985264

Now we combine the two data frames.

fungicide = fungicide_A %>% 
  bind_rows(fungicide_B)
head(fungicide)

  isolate   field   fungicide  dose     growth
1       1 Organic Fungicide A 0e+00 20.2082399
2       1 Organic Fungicide A 1e-05 20.1168279
3       1 Organic Fungicide A 1e-04 19.2479678
4       1 Organic Fungicide A 1e-03 15.8123455
5       1 Organic Fungicide A 1e-02  7.3206757
6       1 Organic Fungicide A 1e-01  0.6985264

There is a way to fit all isolates at once using drm() with curveid = isolate. However, with large datasets this may become slow, and the resulting object can be difficult to handle.

# DO NOT RUN THIS ON VERY LARGE DATASETS WITHOUT CHECKING PERFORMANCE
model = drm(growth~dose, curveid = isolate,  fct = LL.4(), data = fungicide)
ED(model, 50, interval = "delta")

For more practical workflows, I prefer using ec50estimator rather than writing a loop around drm() manually. That is exactly the use case the package was designed for: estimating \(EC_{50}\) isolate by isolate within structured datasets while preserving grouping variables and confidence intervals in a tidy output.

In other words, drc is excellent for understanding and fitting an individual curve, whereas ec50estimator becomes the practical layer when the analysis involves many isolates, fungicides, or production systems.

EC50_est <- estimate_EC50(
  growth ~ dose,
  data = fungicide,
  isolate_col = "isolate",
  strata_col = c("field", "fungicide"),
  interval = "delta",
  fct = drc::LL.4()
) %>% 
  as.data.frame()

head(EC50_est)

  ID   field   fungicide    Estimate   Std..Error       Lower       Upper
1  1 Organic Fungicide A 0.006373092 0.0007047683 0.004935707 0.007810476
2  3 Organic Fungicide A 0.003602134 0.0002580823 0.003075771 0.004128496
3  5 Organic Fungicide A 0.006064882 0.0005081581 0.005028487 0.007101277
4  7 Organic Fungicide A 0.006926374 0.0006407013 0.005619655 0.008233093
5  9 Organic Fungicide A 0.005180691 0.0005403256 0.004078689 0.006282692
6 11 Organic Fungicide A 0.005165555 0.0005897619 0.003962728 0.006368382

Important note: in this post I used the same model, LL.4(), for both fungicides. In real applications, you should verify which model is most appropriate for your data. A useful strategy is to fit candidate models with drm() for a subset of isolates and compare them with mselect(). Once that modeling choice is clear, ec50estimator makes it much easier to apply the workflow consistently across all curves.

Visualizing the estimated EC50 values

Once estimates are available, one useful step is to compare their distributions across strata.

EC50_est %>% 
  ggplot(aes(x = field, y = Estimate, color = field)) +
  geom_jitter(width = 0.1, alpha = 0.6) +
  geom_boxplot(fill = NA, outlier.shape = NA) +
  facet_wrap(~fungicide, scales = "free_x", nrow = 1) +
  theme_half_open(font_size = 12) +
  scale_color_colorblind() +
  scale_y_log10() +
  labs(
    x = "Production system",
    y = expression(EC[50]),
    title = expression("Estimated distributions of " * EC[50] * " across production systems")
  )

Another possibility is to inspect the full distribution more carefully.

EC50_est %>% 
  ggplot(aes(x = Estimate, y = field, fill = field)) +
  geom_density_ridges(alpha = 0.7) +
  facet_wrap(~fungicide, nrow = 2) +
  theme_half_open(font_size = 12) +
  scale_fill_colourblind() +
  scale_x_log10() +
  labs(
    x = expression(EC[50]),
    y = "Production system",
    title = expression("Density view of estimated " * EC[50] * " values")
  )

Warning: `scale_fill_colorblind()` was deprecated in ggthemes 5.2.0.

Picking joint bandwidth of 0.0489

Picking joint bandwidth of 0.0298

Final remarks

Estimating \(EC_{50}\) is an important step in fungicide sensitivity studies, but good inference depends on more than simply fitting one model and reporting one number. It is important to think about:

the biological meaning of the response variable
the appropriateness of the selected dose-response model
the structure of the dataset
how estimates vary across isolates, years, fields, or treatments

The workflow shown here provides a starting point for both teaching and applied analyses in plant pathology.