In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from itertools import product
from scipy.optimize import curve_fit

Подготовка данных

In [2]:
data = pd.read_csv("kinetics_2.csv", sep='\t', index_col=None, header=[0,1])
In [3]:
data.head()
Out[3]:
ABTS, 100 uM ABTS, 208 uM ABTS, 524 uM
time consumption.S time consumption.S time consumption.S
0 0.000 0.00162 0.000 0.00050 0.000 -0.000273
1 0.187 0.02695 0.187 -0.17024 0.250 0.257980
2 0.437 0.45392 0.437 0.35591 0.438 0.637440
3 0.687 0.98718 0.625 0.85406 0.688 0.988040
4 0.875 0.86526 0.875 0.61390 0.938 1.101360
In [4]:
abts_concentrations = [100, 208, 524]
In [5]:
scale_factor = 4 # из уравнения реакции
In [6]:
for conc in abts_concentrations:
    data[(f'ABTS, {conc} uM', 'current_substrate')] = conc - scale_factor * data.loc[:, f'ABTS, {conc} uM']['consumption.S']
In [7]:
data.head()
Out[7]:
ABTS, 100 uM ABTS, 208 uM ABTS, 524 uM ABTS, 100 uM ABTS, 208 uM ABTS, 524 uM
time consumption.S time consumption.S time consumption.S current_substrate current_substrate current_substrate
0 0.000 0.00162 0.000 0.00050 0.000 -0.000273 99.99352 207.99800 524.001092
1 0.187 0.02695 0.187 -0.17024 0.250 0.257980 99.89220 208.68096 522.968080
2 0.437 0.45392 0.437 0.35591 0.438 0.637440 98.18432 206.57636 521.450240
3 0.687 0.98718 0.625 0.85406 0.688 0.988040 96.05128 204.58376 520.047840
4 0.875 0.86526 0.875 0.61390 0.938 1.101360 96.53896 205.54440 519.594560

Задание 1. Уокер-Шмидт.

In [8]:
plt.figure(figsize=(15,5))
for i, conc in enumerate(abts_concentrations):
    x = np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    y = (conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    c = data.loc[:, (f'ABTS, {conc} uM', 'time')]
    #cmap = sns.diverging_palette(240, 10, l=65, center="dark", as_cmap=True)
    cmap = sns.dark_palette("seagreen", reverse=True, as_cmap=True)
    plt.subplot(1,3,i+1)
    plt.scatter(x=x, y=y, c=c, cmap=cmap)
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM')
plt.tight_layout()

Мы видим, что на начальных значениях времени есть линейный участок, причём он начинается с максимума $y$, а кончается на минимуме $x$. Как выбрать оптимальное число первых наблюдений так, чтобы был линейный участок? Будем перебирать и строить графики.

In [9]:
observations = [10, 20, 30, 40, 50]
In [10]:
len(list(product(observations, abts_concentrations)))
Out[10]:
15
In [11]:
plt.figure(figsize=(15, 25))
for i, pair in enumerate(product(observations, abts_concentrations)):
    count, conc = pair
    x = (np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[:count]
    y = ((conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[:count]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')][count]
    #cmap = sns.diverging_palette(240, 10, l=65, center="dark", as_cmap=True)
    plt.subplot(5,3,i+1)
    plt.plot(x, y, 'bo')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM, first {count} observations ({time:.3f} s)')
plt.tight_layout()

Для 100 мкМ лучше всего взять первые 20 наблюдений, для 208 и 524 мкМ подойдёт и 50 наблюдений. Для удобства вычислений будем использовать примерно первые 50 наблюдений для всех концентраций, иначе говоря, первые 10 секунд реакций.

In [12]:
param_list = list()
for conc in abts_concentrations:
    time = data[f'ABTS, {conc} uM'].query('0 < time < 10')['time']
    S = data[f'ABTS, {conc} uM'].query('0 < time < 10')['current_substrate']
    x = np.log(conc / S) / time
    y = (conc - S) / time
    param_list.append(np.polyfit(x, y, deg=1))
In [13]:
print("conc\tKm\tVm")
for conc, i in zip(abts_concentrations, param_list):
    print(f"{conc}\t{-i[0]:.2f}\t{i[1]:.3f}")

conc	Km	Vm
100	-101.32	-0.244
208	-206.17	-0.083
524	-530.27	-0.092

Это очень плохо — отрицательные величины.

In [14]:
plt.figure(figsize=(15,5))
for i, conc in enumerate(abts_concentrations):
    x = np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    y = (conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')]
    plt.subplot(1,3,i+1)
    plt.plot(x, y, 'bo')
    plt.plot(x, param_list[i][0] * x + param_list[i][1], color='red')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM')
plt.tight_layout()

А это значит, что мы выбрали не тот диапазон! Надо смотреть не на первые 10 секунд реакции, а на последующий "линейный" участок. В таком случае, данные при первых секундах реакции можно считать выбросами. Да и точек больше на новом линейном участке.

Найдём начала этих участков.

In [15]:
plt.figure(figsize=(15, 35))
for i, pair in enumerate(product([50, 60, 70, 80, 90, 100, 110], abts_concentrations)):
    count, conc = pair
    x = (np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[count:]
    y = ((conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[count:]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')][count]
    #cmap = sns.diverging_palette(240, 10, l=65, center="dark", as_cmap=True)
    plt.subplot(7,3,i+1)
    plt.plot(x, y, 'bo')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM, first {count} observations ({time:.3f} s)')
plt.tight_layout()

Начиная с 24 секунд пропадают выбросы в левом верхем углу. Вырисовывается линейный участок.

In [22]:
plt.figure(figsize=(15, 25))
for i, pair in enumerate(product([200, 300, 400, 500, 600], abts_concentrations)):
    count, conc = pair
    x = (np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[110:count]
    y = ((conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[110:count]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')][count]
    #cmap = sns.diverging_palette(240, 10, l=65, center="dark", as_cmap=True)
    plt.subplot(5,3,i+1)
    plt.plot(x, y, 'bo')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM, 110 to {count} observations ({time:.3f} s)')
plt.tight_layout()

Для 100 мкМ конец линейного диапазона находится где-то между 200 и 300 наблюдениями (44 секунды и выше). Для 208 мкМ между 400 и 500 (89 и 111 секунд), для 524 мкМ есть излом между 400-ым и 500-ым наблюдениями. Будем перебирать концы для каждой концентрации по отдельности.

In [48]:
plt.figure(figsize=(10, 25))
left = 110
conc = 100
for i, count in enumerate(range(200, 300, 10)):
    x = (np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[left:count]
    y = ((conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[left:count]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')][count]
    #cmap = sns.diverging_palette(240, 10, l=65, center="dark", as_cmap=True)
    plt.subplot(5,2,i+1)
    plt.plot(x, y, 'bo')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM, {left} to {count} observations ({time:.3f} s)')
plt.tight_layout()

Для 100 мкМ область линейной зависимости заканчивается где-то в районе 240-го наблюдения (54 с).

In [49]:
plt.figure(figsize=(10, 25))
left = 110
conc = 208
for i, count in enumerate(range(400, 500, 10)):
    x = (np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[left:count]
    y = ((conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[left:count]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')][count]
    #cmap = sns.diverging_palette(240, 10, l=65, center="dark", as_cmap=True)
    plt.subplot(5,2,i+1)
    plt.plot(x, y, 'bo')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM, {left} to {count} observations ({time:.3f} s)')
plt.tight_layout()

Для 208 мкМ линейный участок заканчивается где-то на 440-ом наблюдении (99 с).

In [50]:
plt.figure(figsize=(10, 25))
left = 110
conc = 524
for i, count in enumerate(range(400, 500, 10)):
    x = (np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[left:count]
    y = ((conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[left:count]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')][count]
    #cmap = sns.diverging_palette(240, 10, l=65, center="dark", as_cmap=True)
    plt.subplot(5,2,i+1)
    plt.plot(x, y, 'bo')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM, {left} to {count} observations ({time:.3f} s)')
plt.tight_layout()

Где-то на 450-ом наблюдении появляется излом (100 с).

Итак, мы получили диапазоны времени для каждой концентрации, где надо смотреть линейные участки.

In [27]:
linear_bounds = {100: [24, 54], 208: [24, 99], 524: [24, 100]}
In [28]:
print("conc\tstart\tstop")
for key, value in linear_bounds.items():
    print(f"{key}\t{value[0]}\t{value[1]}")

conc	start	stop
100	24	54
208	24	99
524	24	100
In [42]:
param_list = list()
for conc in abts_concentrations:
    time = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['time']
    S = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['current_substrate']
    x = np.log(conc / S) / time
    y = (conc - S) / time
    param_list.append(np.polyfit(x, y, deg=1))
In [44]:
print("conc\tKm\tVm")
for conc, i in zip(abts_concentrations, param_list):
    print(f"{conc}\t{-i[0]:.2f}\t{i[1]:.3f}")

conc	Km	Vm
100	21.53	3.219
208	28.20	3.231
524	213.51	4.314

Замечательно, у нас все величиные неотрицательные. Посмотрим на графики.

In [45]:
plt.figure(figsize=(15,5))
for i, conc in enumerate(abts_concentrations):
    x = np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    y = (conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')]
    plt.subplot(1,3,i+1)
    plt.plot(x, y, 'bo')
    plt.plot(x, param_list[i][0] * x + param_list[i][1], color='red')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM')
plt.tight_layout()

Видим, что для 100 и 208 мкМ данные по $K_M$ примерно согласуются между собой и хорошо фитируют линейный участок. Для 524 мкМ это сделано плохо, $K_M$ там огромная. Заново найдём линейный участок для 524 мкМ. Мы знаем, что излом начинается на 450-ом наблюдении, от него и будем искать конец линейного участка.

In [52]:
len(range(450, 1300, 50))
Out[52]:
17
In [53]:
plt.figure(figsize=(20, 25))
left = 450
conc = 524
for i, count in enumerate(range(450, 1300, 50)):
    x = (np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[left:count]
    y = ((conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')])[left:count]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')][count]
    #cmap = sns.diverging_palette(240, 10, l=65, center="dark", as_cmap=True)
    plt.subplot(5,4,i+1)
    plt.plot(x, y, 'bo')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM, {left} to {count} observations ({time:.3f} s)')
plt.tight_layout()

Так как при 524 мкМ есть два участка линейности (100 — 156 секунд и 156 — 245 секунд), то неясно, какой из них выбрать.

Пусть выберем 100 — 156 секунд.

In [67]:
linear_bounds = {100: [24, 54], 208: [24, 99], 524: [100, 156]}
print("conc\tstart\tstop")
for key, value in linear_bounds.items():
    print(f"{key}\t{value[0]}\t{value[1]}")

conc	start	stop
100	24	54
208	24	99
524	100	156
In [68]:
param_list = list()
for conc in abts_concentrations:
    time = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['time']
    S = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['current_substrate']
    x = np.log(conc / S) / time
    y = (conc - S) / time
    param_list.append(np.polyfit(x, y, deg=1))
print("conc\tKm\tVm")
for conc, i in zip(abts_concentrations, param_list):
    print(f"{conc}\t{-i[0]:.2f}\t{i[1]:.3f}")

conc	Km	Vm
100	21.53	3.219
208	28.20	3.231
524	131.49	3.705
In [69]:
plt.figure(figsize=(15,5))
for i, conc in enumerate(abts_concentrations):
    x = np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    y = (conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')]
    plt.subplot(1,3,i+1)
    plt.plot(x, y, 'bo')
    plt.plot(x, param_list[i][0] * x + param_list[i][1], color='red')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM')
plt.tight_layout()

Всё равно получаются несогласующиеся результаты по $K_M$, $V_m$ ещё кое-как сходится.

Выберем отрезок 156 — 245 секунд.

In [71]:
linear_bounds = {100: [24, 54], 208: [24, 99], 524: [156, 245]}
print("conc\tstart\tstop")
for key, value in linear_bounds.items():
    print(f"{key}\t{value[0]}\t{value[1]}")

conc	start	stop
100	24	54
208	24	99
524	156	245
In [72]:
param_list = list()
for conc in abts_concentrations:
    time = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['time']
    S = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['current_substrate']
    x = np.log(conc / S) / time
    y = (conc - S) / time
    param_list.append(np.polyfit(x, y, deg=1))
print("conc\tKm\tVm")
for conc, i in zip(abts_concentrations, param_list):
    print(f"{conc}\t{-i[0]:.2f}\t{i[1]:.3f}")

conc	Km	Vm
100	21.53	3.219
208	28.20	3.231
524	40.61	2.861
In [73]:
plt.figure(figsize=(15,5))
for i, conc in enumerate(abts_concentrations):
    x = np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    y = (conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')]
    plt.subplot(1,3,i+1)
    plt.plot(x, y, 'bo')
    plt.plot(x, param_list[i][0] * x + param_list[i][1], color='red')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM')
plt.tight_layout()

$K_M$ уже не очень сильно отличается от других, но $V_m$ как-то просела.

Что, если выбрать сразу оба диапазона?

In [74]:
linear_bounds = {100: [24, 54], 208: [24, 99], 524: [100, 245]}
print("conc\tstart\tstop")
for key, value in linear_bounds.items():
    print(f"{key}\t{value[0]}\t{value[1]}")

conc	start	stop
100	24	54
208	24	99
524	100	245
In [75]:
param_list = list()
for conc in abts_concentrations:
    time = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['time']
    S = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['current_substrate']
    x = np.log(conc / S) / time
    y = (conc - S) / time
    param_list.append(np.polyfit(x, y, deg=1))
print("conc\tKm\tVm")
for conc, i in zip(abts_concentrations, param_list):
    print(f"{conc}\t{-i[0]:.2f}\t{i[1]:.3f}")

conc	Km	Vm
100	21.53	3.219
208	28.20	3.231
524	53.02	3.036
In [76]:
plt.figure(figsize=(15,5))
for i, conc in enumerate(abts_concentrations):
    x = np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    y = (conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')]
    plt.subplot(1,3,i+1)
    plt.plot(x, y, 'bo')
    plt.plot(x, param_list[i][0] * x + param_list[i][1], color='red')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM')
plt.tight_layout()

$K_M$ стала сильнее отличаться от остальных, $V_m$ почти попала.

Всё же надо идти ещё дальше по времени реакции, начинать не со 156 секунд, а позже, например, со 200 секунд.

In [87]:
linear_bounds = {100: [24, 54], 208: [24, 99], 524: [200, 245]}
print("conc\tstart\tstop")
for key, value in linear_bounds.items():
    print(f"{key}\t{value[0]}\t{value[1]}")

conc	start	stop
100	24	54
208	24	99
524	200	245
In [88]:
param_list = list()
for conc in abts_concentrations:
    time = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['time']
    S = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['current_substrate']
    x = np.log(conc / S) / time
    y = (conc - S) / time
    param_list.append(np.polyfit(x, y, deg=1))
print("conc\tKm\tVm")
for conc, i in zip(abts_concentrations, param_list):
    print(f"{conc}\t{-i[0]:.2f}\t{i[1]:.3f}")

conc	Km	Vm
100	21.53	3.219
208	28.20	3.231
524	32.86	2.741
In [89]:
plt.figure(figsize=(15,5))
for i, conc in enumerate(abts_concentrations):
    x = np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    y = (conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')]
    plt.subplot(1,3,i+1)
    plt.plot(x, y, 'bo')
    plt.plot(x, param_list[i][0] * x + param_list[i][1], color='red')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM')
plt.tight_layout()

Получили $K_M$ более-менее приличный, но по $V_m$ упали. Оставим такие результаты, видимо, данные не самые хорошие.

Задание 2. Нелинейная оптимизация.

In [90]:
def integral_func(x, km, vm):
    return km / vm * (np.log(x[:, 0]) - np.log(x[:, 1])) + 1 / vm * (x[:, 0] - x[:, 1])
In [93]:
param_list = list()
for conc in abts_concentrations:
    S = data[f'ABTS, {conc} uM']['current_substrate'].dropna()
    x = np.array(list(product([conc], S)))
    y = data[f'ABTS, {conc} uM']['time'].dropna()
    params, cov = curve_fit(integral_func, x, y)
    param_list.append(params)
In [95]:
print("conc\tKm\tVm")
for conc, i in zip(abts_concentrations, param_list):
    print(f"{conc}\t{i[0]:.2f}\t{i[1]:.3f}")

conc	Km	Vm
100	27.50	3.463
208	37.47	3.451
524	55.27	3.070

Весьма странно, что $K_M$ увеличивается с ростом концентрации, но пусть так. $V_m$ примерно одинакова. Данные для 524 мкМ показывают, что для линейной аппроксимации надо было выбрать диапазон 100 — 245 секунд, тогда результаты примерно сходятся.

In [103]:
plt.figure(figsize=(15,5))
for i, conc in enumerate(abts_concentrations):
    S = data[f'ABTS, {conc} uM']['current_substrate'].dropna()
    x = np.array(list(product([conc], S)))
    y = data[f'ABTS, {conc} uM']['time'].dropna()
    plt.subplot(1,3,i+1)
    plt.plot(S, y, 'bo')
    plt.plot(S, integral_func(x, *param_list[i]), color='red')
    plt.xlabel("current S")
    plt.ylabel("Time, sec")
    plt.title(f'ABTS, {conc} uM')
plt.tight_layout()

Кривые зафитились прекрасно.

In [102]:
plt.figure(figsize=(15,5))
for i, conc in enumerate(abts_concentrations):
    x = np.log(conc / data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    y = (conc - data.loc[:, (f'ABTS, {conc} uM', 'current_substrate')]) / data.loc[:, (f'ABTS, {conc} uM', 'time')]
    time = data.loc[:, (f'ABTS, {conc} uM', 'time')]
    plt.subplot(1,3,i+1)
    plt.plot(x, y, 'bo')
    plt.plot(x, - param_list[i][0] * x + param_list[i][1], color='red')
    plt.xlabel(r"$\frac{ln(\frac{S_0}{S})}{t}$", fontsize='x-large')
    plt.ylabel(r"$\frac{(S_0 - S)}{t}$", fontsize='x-large')
    plt.title(f'ABTS, {conc} uM')
plt.tight_layout()

Если посмотреть на график для первого задания с $K_M$ и $V_m$ из второго задания, то видно, что результаты более-менее согласуются между двумя методами. Данные для 524 мкМ не обладают той же линейностью, что и для других концентраций. У 100 и 208 мкМ есть лишь один линейный участок между началом и концом наблюдений, а у 524 мкМ там кривая.

Сравнение параметров из разных аппроксимаций.

In [104]:
param_linear = list()
linear_bounds = {100: [24, 54], 208: [24, 99], 524: [100, 245]}
param_list = list()
for conc in abts_concentrations:
    time = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['time']
    S = data[f'ABTS, {conc} uM'].query(f'{linear_bounds[conc][0]} < time < {linear_bounds[conc][1]}')['current_substrate']
    x = np.log(conc / S) / time
    y = (conc - S) / time
    param_linear.append(np.polyfit(x, y, deg=1))
print("conc\tKm\tVm")
for conc, i in zip(abts_concentrations, param_linear):
    print(f"{conc}\t{-i[0]:.2f}\t{i[1]:.3f}")

conc	Km	Vm
100	21.53	3.219
208	28.20	3.231
524	53.02	3.036
In [106]:
param_nonlinear = list()
for conc in abts_concentrations:
    S = data[f'ABTS, {conc} uM']['current_substrate'].dropna()
    x = np.array(list(product([conc], S)))
    y = data[f'ABTS, {conc} uM']['time'].dropna()
    params, cov = curve_fit(integral_func, x, y)
    param_nonlinear.append(params)
print("conc\tKm\tVm")
for conc, i in zip(abts_concentrations, param_nonlinear):
    print(f"{conc}\t{i[0]:.2f}\t{i[1]:.3f}")

conc	Km	Vm
100	27.50	3.463
208	37.47	3.451
524	55.27	3.070
In [117]:
linear_array = np.array(param_linear)
linear_array[:, 0] *= -1
nonlinear_array = np.array(param_nonlinear)
params_df = pd.DataFrame(np.concatenate([linear_array, nonlinear_array], axis=0)).rename(columns={0: 'Km', 1: 'Vm'})
In [119]:
params_df['method'] = ['linear' for _ in range(3)] + ['nonlinear' for _ in range(3)]
params_df['concentration'] = abts_concentrations + abts_concentrations
In [120]:
params_df
Out[120]:
Km Vm method concentration
0 21.528077 3.219087 linear 100
1 28.196821 3.231493 linear 208
2 53.022352 3.036234 linear 524
3 27.498343 3.462998 nonlinear 100
4 37.468611 3.450596 nonlinear 208
5 55.265330 3.069645 nonlinear 524
In [127]:
plt.figure(figsize=(10, 5))
plt.subplot(1,2,1)
sns.scatterplot(data=params_df, x='concentration', y='Km', hue='method')
plt.subplot(1,2,2)
sns.scatterplot(data=params_df, x='concentration', y='Vm', hue='method')
plt.tight_layout()

Обе аппроксимации выявляют зависимость $K_M$ и $V_m$ от исходной концентрации субстрата, отличаются лишь некоторым сдвигом в данных (при 528 мкМ он минимален). Различие обусловлено тем, что линейную аппроксимацию мы обучаем на куске данных, а нелинейную -- на всём объёме данных, так что там получаем дополнительную информацию (или шум).

Метод нелинейной аппроксимации функцией $t = \frac{K_M}{V_m} \ast ln \frac{S_0}{S} + \frac{1}{V_m} \ast (S_0 - S)$ более точно зафитировался на наших данных, чем линейная аппроксимация в координатах Уокера-Шмидта. Это связано с тем, что в координатах Уокера-Шмидта линеаризация оказалась не самой хорошей.