Konzentrations-Wirkungsbeziehungen am Beispiel von Antibiotika
1 Motivation
1.1 Inhalt der Übung
In dieser Übung erstellen wir Konzentrations-Wirkungsbeziehungen. Als toxische Substanzen dienen Antibiotika und als Testorganismen unterschiedliche Stämme von Escherichia coli. Aufgrund der außerordentlich großen genetischen Variabilität von Bakterien, rechnen wir mit Unterschieden zwischen den getesteten Isolaten. Als konkretes Kriterium zum Messen der Fitness (a.k.a. ‘end point’) verwenden wir die Wachstumsrate.
1.2 Praktische Bedeutung
Konzentrations-Wirkungsbeziehungen sind allgemein nützlich, um die Wirkung von Substanzen knapp und eindeutig mathematisch zu beschreiben. Man kann sie bspw. verwenden, um die Folgen einer Exposition gegenüber toxischen Substanzen abzuschätzen (mittels mathematischer Modellierung) oder notwendige Dosen für Therapien/Bekämpfungsmaßnahmen festzulegen.
2 Theoretische Grundlagen
2.1 Definition der Wachtsumsrate
In Phasen ohne Ressourcenlimitation wächst eine Kultur von Mikroorganismen exponentiell in Folge fortwährender Teilung. Die Entwicklung der Dichte der Kultur (\(N\), Zellen / mL) wird dann durch eine gewöhnliche Differentialgleichung beschrieben (Equation 1). Darin hat die Wachstumsrate \(\mu\) üblicherweise die Einheit 1/h.
\[ \frac{d}{dt} N = \mu \cdot N \tag{1}\]
Die Lösung von Equation 1 bei bekanntem Anfangswert \(N(t_0)\) ist eine exponentielle Funktion der Zeit (Equation 2).
\[ N(t_0 + \Delta t) = N(t_0) \cdot exp(\mu \cdot \Delta t) \tag{2}\]
Durch Umstellen von Equation 2 erhält man eine Vorschrift zur Schätzung der Wachstumsrate (Equation 3).
\[ \mu = \frac{ln\left( \frac{N(t_0 + \Delta t)}{N(t_0)} \right)}{\Delta t} \tag{3}\]
Um \(\mu\) zu bestimmen, werden also mindestens zwei Schätzungen der Kulturdichte zu unterschiedlichen Zeitpunkten innerhalb der exponentiellen Wachstumsphase benötigt. Ob sich eine Kultur in der exponentiellen Phase befindet oder nicht, lässt sich aus einem Plot der Dynamik mit logarithmierter y-Achse unmittelbar erkennen (Figure 1). Die exponentielle Phase ist dort als linearer Bereich erkennbar. Dies folgt unmittelbar aus Equation 2, denn wenn diese logarithmiert wird, ist die rechte Seite eine lineare Funktion der Zeit (Equation 4).
\[ ln \left( N(t_0 + \Delta t) \right) = ln \left( N(t_0) \right) + \mu \cdot \Delta t \tag{4}\]
2.2 Messgröße
Um Wachstumsraten zu messen, müssen wir die Entwicklung der Dichte der Kultur über die Zeit aufzeichnen (Populationsdynamik). Wir verwenden dabei die Optische Dichte (a.k.a. Extinktion) als Proxy für die Zelldichte. Über einen weiten Bereich besteht zwischen optischer Dichte und Zelldichte ein linearer Zusammenhang gemäß dem Lambert-Beer’schen Gesetz. Üblichweise wird die optische Dichte bei 600 nm gemessen (oranges Licht) und kurz als OD600 bezeichnet.
Aus Equation 3 ist ersichtlich, dass wir zur Bestimmung der Wachstumsrate lediglich einen Quotienten zweier Dichtemessungen bzw. den Anstieg der gerade im rechten Teil von Figure 1 auswerten. Deshalb ist, im Gegensatz zu sonst üblichen Anwendungen der Photometrie, keine Eichkurve nötig, sofern nur \(\mu\) bestimmt werden soll.
2.3 Strukturierung der Daten
Damit die spätere Auswertung leicht fällt, müssen die Daten wohl organisiert gespeichert werden. Im konreten Fall haben wir lediglich mit einer Messgröße zu tun (OD600) aber fünf weitere Variablen zu berücksichtigen:
- die Kennung des Isolats, sofern mehrere getestet werden
- das verwendete Antibiotikum, sofern mehrere getestet werden
- die Konzentration des Antibiotikums, da wir einen Gradienten untersuchen
- die Kennung des Replikats, sofern physische Parallelen untersucht werden
- den Zeitpunkt der Messung, z.B. in Stunden nach Beginn des Experiments
In der Praxis bietet es sich an, die Daten unmittelbar so abzulegen, wie sie in einer Datenbank gespeichert würden, d.h. als eine Tabelle mit 6 Spalten und (sehr) vielen Zeilen.
isolate | antibiotic | conc | replicate | hour | OD600
--------+------------+------+-----------+------+------
| | | | |
| | | | |
| | | | |Die Benutzung eines Tabellenkalkulationsprogramms bietet sich an, da man Blöcke sich wiederholender Information bequem ausfüllen kann. Aber Vorsicht: Dabei passieren auch leicht Fehler. Ein leserliches Protokoll auf Papier sollte auf jeden Fall geführt, bis zum Schluss aufbewahrt, und als backup fotographiert werden.
3 Datenanalyse
3.1 Überblick
Unsere Datenanalyse umfasst eine Reihe aufeinander aufbauender Schritte. Wie üblich, geht es los mit:
- Importieren der Daten
- Vorverarbeitung / Bereinigen problematischer Werte
- Visualisieren
Anschließen werden wir:
- … für jede Kombination von Antibiotikum, Konzentration, Isolat, und Replikat eine Wachstumsrate schätzen,
- … für jede Kombination von Antibiotikum und Isolat ein statistisches Modell fitten, welches die Abhängigkeit des Wachstums von der Konzentration beschreibt.
In den Schritten 4 und 5 reduziert sich jeweils die Komplexität der Daten, repräsentiert durch die Anzahl der benötigten Tabellenspalten, wie folgt:
Vor Schritt 4 | Nach Schritt 4 | Nach Schritt 5
-----------------+------------------------+------------------
Isolat | Isolat | Isolat
Antibiotikum | Antibiotikum | Antibiotikum
Konzentration | Konzentration --+--> | K.-W.-Beziehung
Replikat | Replikat --+ |
Stunde --+--> | Wachstumsrate --+ |
OD600 --+ |Die Datenanalyse lässt sich mittels R in unterschiedlichen Programmierstilen implementieren. Ich benutze hier die traditionelle Variante. Zentrales Element ist die unscheinbare Funktion by. Diese teilt die in einem data.frame abgelegten Daten in Scheiben (Gruppen von Zeilen) und wendet auf jede der Scheiben eine Funktion an.
3.2 Daten einlesen und vorbereiten
Wir lesen die Daten aus einer Excel-Datei im “long-format”, so wie oben vorgeschlagen. Dateiname und Arbeitsverzeichnis sind ggf. anzupassen. Der verwendete Testdatensatz enthält keine physischen Replikate, weshalb die entsprechende Kennung durchgehend “1” lautet.
rm(list=ls())
library("readxl")
#x <- as.data.frame(read_excel("referenceDataset.xlsx"))
x <- as.data.frame(read_excel("students2026Dataset.xlsx"))
x[x[,"antibiotic"] == "amp", "antibiotic"] <- "Ampicillin"
x[x[,"antibiotic"] == "tet", "antibiotic"] <- "Tetracycline"Und so sieht der Anfang der Tabelle aus:
isolate antibiotic conc replicate hour OD600
1 A Ampicillin 12.5000 HK 0 0.001
2 A Ampicillin 6.2500 HK 0 0.001
3 A Ampicillin 3.1250 HK 0 0.001
4 A Ampicillin 1.5625 HK 0 0.001
5 A Ampicillin 0.0000 HK 0 0.001
6 A Ampicillin 12.5000 HK 1 0.008In den Daten gibt es einige leicht negative Werte, d.h. die Proben waren scheinbar weniger trüb als die Blindprobe ohne Bakterien (blank). Absolut sind die Werte allerdings sehr gering und im Bereich der Auflösung des Messgeräts. Denoch stören sie, da wir die Daten später logarithmieren müssen. Als pragmatische Lösung heben wir alle Wert eines konkreten Datensatzes so an, dass der kleinste Messwert bei 0.001 liegt.
shift <- function(x) {
if (any(x[,"OD600"] <= 0)) {
x[,"OD600"] <- x[,"OD600"] + (0.001 - min(x[,"OD600"]))
}
x
}
x <- do.call(rbind, by(x, x[,c("antibiotic","isolate","replicate")], shift))3.3 Visualisierung
Zuerst plotten wir die Beobachtungen über die Zeit, getrennt nach Isolat und Antibiotikum (Figure 2). Die Konzentrationen lassen sich durch Farben kodieren. Durch die logarithmische y-Achse, kann man sofort die Phase des exponentiellen Wachstums (und ggf. Abweichungen davon) erkennen.
plot_dynamics <- function(x) {
clr <- colorRampPalette(c("steelblue4","skyblue","khaki","red"))(
length(unique(x[,"conc"])))
plot(x[,"hour"], x[,"OD600"], log="y", bty="n", yaxt="n",
xlim=c(0, max(x[,"hour"])), xlab="Hour", ylab="",
col=clr[match(x[,"conc"], sort(unique(x[,"conc"])))], pch=20)
axis(side=2, las=2)
mtext(side=2, line=4, cex=par("cex"),
expression(paste("OD",{}[600]," (log scaled)")))
legend("topleft", bty="n", ncol=2, pch=20, col=clr,
legend=sort(round(unique(x[,"conc"]),2)))
mtext(side=3, at=0, adj=0, cex=par("cex"),
paste0("Isolate '",unique(x[,"isolate"]),"' with ",
unique(x[,"antibiotic"])))
invisible(NULL)
}
par(mfcol=c(length(unique(x[,"isolate"])), length(unique(x[,"antibiotic"]))))
omar <- par("mar")
par(mar=c(4.5, 5.5, 2, 1))
ignored <- by(x, x[c("isolate","antibiotic")], plot_dynamics)
par(mar=omar)
par(mfcol=c(1,1))
3.4 Schätzung der (maximalen) Wachstumsraten
Als Indikator für die Fitness wollen wir die maximalen Wachstumsraten (\(\mu_{max}\)) berechnen. Dafür suchen wir die Phase des steilsten Anstiegs der logarithmierten OD-Wert und fitten dort ein lineares Modell (vgl. Equation 4). Bei der Auswahl der Werte in der Phase des steilsten Anstiegs gibt es einiges zu beachten:
Wir müssen uns im zuverlässigen Messbereich des Photometers befinden. Das ist ganz zu Beginn des Experiments nicht unbedingt der Fall. Wenn das Photometer z.B. drei Nachkommastellen ausgibt, sind Werte wie 0.001 oder 0.003 mit sehr großen Unsicherheiten (im Sinne von relativen Fehlern) behaftet.
Wir sollten die ausgewertete Zeitspanne so lang wählen, dass wir mindestens 3 Zeitpunkte (= 2 Intervalle) verarbeiten. Andererseits dürfen wir nur Zeitpunkte aus der exponentiellen Phase wählen, und da bleibt uns aufgrund des schnellen Wachstums in der Praxis wenig Auswahl (vgl. Figure 1).
Für das fitten von Wachstumsraten gibt es unterschiedliche Ansätze und auch mehrere fertige R-Pakete, wie z.B. dieses. Wir benutzen hier einen einfachen, sehr robusten Ansatz, der sich leicht implementieren und nachvollziehen lässt (siehe Funktion mu_max).
mu_max <- function(
x, # table with growth curve data (single case)
time="time", # column with time information
dens="OD600", # column with culture density information
intervals=2, # number of time intervals to consider
min_time=0.5 # drop times earlier than this
) {
stopifnot(all(c(time, dens) %in% colnames(x)))
x <- x[x[,time] >= min_time,]
stopifnot(intervals < nrow(x))
stopifnot(length(unique(x[,time])) == nrow(x))
x <- x[order(x[,time]),]
rng.try <- 1:(intervals+1)
mu <- -Inf
while (max(rng.try) < nrow(x)) {
fit <- lm(log(y) ~ time,
data=data.frame(time=x[rng.try,time], y=x[rng.try,dens]))
mu <- max(mu, coef(fit)["time"])
rng.try <- rng.try + 1
}
mu
}
splitCols <- c("isolate","antibiotic", # columns to split cases
"conc","replicate")
x <- by(x, x[,splitCols], mu_max, # process all cases
time="hour", dens="OD600")
x <- array2DF(x, responseName="mu_max") # reformat output
x[,"conc"] <- as.numeric(x[,"conc"])
x <- x[is.finite(x[,"mu_max"]),] # drop non-existing casesUnd so sieht der Anfang der Ergebnistabelle aus. Die Werte von \(\mu_{max}\) sind in der Einheit 1/Stunde. Für die leichtere Interpretation ist auch die Verdopplungszeit \(\tau\) in Minuten mit angegeben. Bekannt ist, dass \(\tau\) für E. coli unter optimalen Bedingungen bei etwa 18 Minuten liegt.
isolate antibiotic conc replicate mu_max tau
4 D Ampicillin 0.000 A. B. M. 1.3849906 30
8 D Tetracycline 0.000 A. B. M. 1.3442071 31
16 D Tetracycline 0.195 A. B. M. 0.7702225 54
24 D Tetracycline 0.390 A. B. M. 0.3209269 130
32 D Tetracycline 0.780 A. B. M. 0.2027326 205
40 D Tetracycline 1.560 A. B. M. 0.1256572 3313.5 Berechnen eines dimensionlosen Fitness-Index
Es ist praktisch, \(\mu_{max}\) zunächst auf den Wert des Kontrollexperiments ohne Antibiotikum zu normieren. So erhält man einen dimensionslosen Indikator für die Fitness, dessen Wertebereich zwischen 0 (volle Hemmung) und 1 (keine Hemmung) liegt.
# find controls
ctrl <- x[(x[,"conc"] == 0), names(x) != "conc"]
# assign controls to records
x <- merge(x, ctrl, by=c("isolate", "antibiotic", "replicate"),
suffixes=c("","_ctrl"))
# normalize
x <- cbind(x, fitness= x[,"mu_max"] / x[,"mu_max_ctrl"])
print(head(x)) isolate antibiotic replicate conc mu_max tau mu_max_ctrl tau_ctrl fitness
1 A Ampicillin AH 0.0000 1.3451923 31 1.345192 31 1.0000000
2 A Ampicillin AH 6.2500 1.2264687 34 1.345192 31 0.9117423
3 A Ampicillin AH 12.5000 0.6479755 64 1.345192 31 0.4816973
4 A Ampicillin AH 1.5625 1.3926003 30 1.345192 31 1.0352425
5 A Ampicillin AH 3.1250 1.3517199 31 1.345192 31 1.0048526
6 A Ampicillin ER 0.0000 1.8863805 22 1.886380 22 1.00000003.6 Fitten von Konzentrations-Wirkungsbeziehungen
3.6.1 Warum Modelle anpassen?
Die Fitness setzen wir nun in Beziehung zur Konzentration, natürlich separat für jedes Isolat und jedes Antibiotikum. An den Zusammenhang zwischen Konzentration und Fitness (summarisch über alle Replikate) können wir ein Modell anpassen, um die Konzentrations-Wirkungsbeziehung mathematisch zu beschreiben. Das macht Sinn, weil …
man die Modelle zur Interpolation und Extrapolation benutzen kann.
man die Parameter der Modelle einfach vergleichen kann, was mit Rohdaten schwierig ist.
Ein üblicher Parameter für derartige Vergleiche ist z.B. die EC50, als Konzentation, welche den halbmaximalen Effekt zeigt. Für das vorliegende Beispiel ist jedoch die minimale Hemmkonzentration eher geeignet und üblich.
3.6.2 Modellanpassung braucht Erfahrung
Das Anpassen der Modelle ist in der Praxis oft nicht trivial. Das hat 2 Gründe:
Wir müssen ein Modell auswählen, was zu den Daten und der Theorie passt.
Wir müssen die Parameter des Modells schätzen. Das ist nur bei linearen Modellen trivial, erfordert ansonsten aber eine Optimierung, mit allen Tücken (z.B. fehlende Konvergenz oder lokale Minima).
3.6.3 Auswahl des Modells
Um ein passendes Modell auszuwählen, könnten wir natürlich einfach kopieren, was andere schon gemacht haben oder dem Vorschlag der KI folgen. Ich empfehle aber, selbst nach- oder mindestens mitzudenken. Klar ist, dass das Modell bei einer Konzentration von Null eine Fitness von Eins zurückgeben muss. Insbesondere sollte man also die Frage klären, wie sich ein plausibles Modell bei sehr hohen Konzentrationen verhalten müsste:
Erwarten wir eine (langsame) Konvergenz der Fitness gegen Null? Falls ja, kommt ein sigmoidales Modell (S-Kurve) oder eine stetig fallende Potenzfunktion in Frage.
Oder sollte die Fitness bei einem Schwellenwert tatsächlich Null erreichen und bei weiter steigenden Konzentrationen bei Null verharren? Dann kommen selbst lineare Modelle in Betracht.
Sofern die Abbildung hormetischer Effekte gewünscht wäre (hier nicht der Fall), müssten Funktionen verwendet werden, die ein Maximum abbilden können.
3.6.4 Praktische Anpassung
Für den gegebenen Datensatz wird pragmatisch ein lineares Modell (Equation 5) und ein einfaches nicht-lineares Modell (Equation 6) angepasst. Beide Modelle erfüllen die Bedingung, dass die Fitness (linke Seite) bei einer Konzentration von Null den Wert Eins annimmt und bei Erreichen der MIC den Wert Null.
\[ \frac{\mu_{max}}{\mu_{max}(0)} = 1 - \frac{conc}{MIC} \tag{5}\]
\[ \frac{\mu_{max}}{\mu_{max}(0)} = 1 - \left( \frac{conc}{MIC} \right)^b \tag{6}\]
Wir definieren zunächst eine Funktion, welche die Anpassungen für eine konkrete “Scheibe” der Daten, d.h. jede Kombination von Isolat und Antibiotikum, vornimmt. Bei der Anpassung des linearen Modells ist zu beachten, dass wir hier einen y-Achsenabschnitt (intercept) von 1 fest vorgeben wollen, was etwas umständlich ist. Wir subtrahieren hier zunächst die 1 von den Daten, passen ein Modell mit dem Intercept Null an, und addieren die 1 wieder.
Für die Anpassung des nichtlinearen Modells könnte man auch die vorgefertigte Funktion nls verwenden. Ich nutze hier, der Robustheit, die low-level Variante mittels optim.
conc_resp <- function(
x, # table of fitness at varying concentrations
min_fitness=0.1 # exclude from fitting if fitness is < this
) {
x <- x[order(x[,"conc"]),]
plot(x[,"conc"], x[,"fitness"], ylim=c(0, 1), bty="L", pch=20,
xlab="mg/L", ylab=expression(paste("µ",{}[max]," / µ",{}[max],"(0)")))
abline(h=seq(0, 1, 0.2), lty=3, col="grey")
legend("bottomleft", bty="n", legend=paste0("'",
unique(x[,"isolate"]),"' with ",unique(x[,"antibiotic"])))
# fit model
x <- x[x[,"fitness"] >= min_fitness,]
out <- ""
## linear without intercept
fit <- lm(I(fitness-1) ~ conc + 0, data=x)
lines(x[,"conc"], predict(fit, newdata=x) + 1, col="red")
mic_linear <- -1 / coef(fit)[["conc"]]
## non-linear
fitness <- function(p, conc) { with(as.list(p), (1 - (conc/mic)^b)) }
ssr <- function(p, x) { sum((fitness(p, x[,"conc"]) - x[,"fitness"])^2) }
fit <- optim(c(mic=mic_linear, b=1), ssr, x=x)
if (fit$convergence == 0) {
test.conc <- seq(0, max(x[,"conc"]), length.out=20)
lines(test.conc, fitness(fit$par, test.conc), col="blue")
mic_nonlinear <- fit$par["mic"]
} else {
mic_nonlinear <- NA
}
mtext(side=3, cex=par("cex"), paste0("MIC: linear = ",
signif(mic_linear, 2),", non-linear = ", signif(mic_nonlinear, 2)))
data.frame(
using_replicates = length(unique(x[,"replicate"])) > 1,
isolate = unique(x[,"isolate"]),
antibiotic = unique(x[,"antibiotic"]),
MIC_linear = signif(mic_linear,2),
MIC_nonlinear = signif(mic_nonlinear,2),
row.names=NULL
)
}Nachfolgend wenden wir die Funktion auf alle Datenscheiben an (Figure 3). Sofern der Datensatz für jede Kombination von Isolat und Antibiotikum mehrere Replikate enthält, werden diese sämtlich genutzt.
par(mfcol=c(length(unique(x[,"isolate"])), length(unique(x[,"antibiotic"]))))
omar <- par("mar")
par(mar=c(4.5, 5.5, 2, 1))
splitCols <- c("isolate","antibiotic") # columns to split cases
out <- by(x, x[c("isolate","antibiotic")], conc_resp) # process all cases
par(mar=omar)
par(mfcol=c(1,1))
In der Praxis enthalten Datensätze häufig Ausreißer. Die Gründe sind vielfältig und reichen von Pipettierfehlern, über Fingerabdrücke auf der Küvette bis zu Schreibfehlern in der Aufzeichnung. Sofern mehrere Replikate vorliegen, können wir eine erneute Anpassung von Modellen versuchen, nachdem Extremwerte entfernt wurden. Ein einfaches Mittel besteht darin, lediglich Medianwerte anstelle der Information einzelner Replikate zu verwenden.
# replace individual fitness estimates by the median
fn <- function(x) {
cbind(unique(x[,c("isolate","antibiotic","conc")]),
replicate="median", fitness=median(x[,"fitness"]))
}
x <- do.call(rbind, by(x, x[,c("isolate","antibiotic","conc")], fn))Die Anpassung der Modelle erfolgt dann wie zuvor.
par(mfcol=c(length(unique(x[,"isolate"])), length(unique(x[,"antibiotic"]))))
omar <- par("mar")
par(mar=c(4.5, 5.5, 2, 1))
splitCols <- c("isolate","antibiotic") # columns to split cases
out <- by(x, x[c("isolate","antibiotic")], conc_resp) # process all cases
par(mar=omar)
par(mfcol=c(1,1))
Letztlich können wir die Schätzwerte noch übersichtlich als Tabelle ausgeben lassen.
print(do.call(rbind, out)) using_replicates isolate antibiotic MIC_linear MIC_nonlinear
1 FALSE A Ampicillin 35.00 20
2 FALSE B Ampicillin 28.00 33
3 FALSE C Ampicillin 36.00 44
4 FALSE D Ampicillin -6800.00 -6800
5 FALSE A Tetracycline 1.50 2
6 FALSE B Tetracycline 58.00 58
7 FALSE C Tetracycline 3.20 660
8 FALSE D Tetracycline 0.77 14 Interpretation
Da die Interpretation weitgehend selbständig erfolgen soll, werden hier nur wenige Tips gegeben.
4.1 Plausibilität
Im Fall von Ampicillin ist in Figure 2 zu erkennen, dass der Effekt keine sofortige Wachstumshemmung ist, sondern erst verzögert einsetzt und eher einen bakteriziden statt einen bakteriostatischen Charakter hat (sofern keine andere Ursache vorliegt). Im Fall bakterizider Antibiotika (bei höheren Dosen) würde man statt Wachstum eher die Mortalität auswerten.
4.2 Einordnung der geschätzten MICs
Vergleichswerte findet man in der EUCAST Datenbank für eine große Zahl von Spezies und Antibiotika. Aufgeführt ist dort für die Anzahl der getesteten Isolate, deren MIC der im Tabellenkopf angegebenen Konzentrationsstufe entsprach.
Für E. coli und die Antibiotika Ampicillin und Tetracycline ergeben sich jeweils bimodale Verteilungen (Figure 5). Dabei repräsentiert die linke Komponente der Verteilung sensible Isolate (geringe MIC, in gelb), die rechte Komponente repräsentiert resistente Isolate bzw. solche mit erhöhter Toleranz (in blau). Der Übergang zwischen den beiden Komponenten wird von EUCAST als Schwellenwerte namens ECOFF angegeben. Er beträgt hier in beiden Fällen 8 mg/L.
Beim Vergleich der geschätzten MICs (Figure 4) mit der statistischen Verteilung (Figure 5) muss bedacht werden, dass beiden Datensätzen unterschiedliche Methoden zugrunde liegen. Eine Einordnung der Sensibilität bzw. Resistenz der getesten Isolate ist aber möglich.