import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import curve_fit
from itertools import product
Загрузка данных¶
Фермент: простагландин-H-синтаза (PGHS) из везикулярных желёз барана (Ovis aries).
Варьируемый субстрат: арахидоновая кислота.
Ингибитор: напроксен.
data = pd.read_csv("2019_практикум-3_302_Мыларщиков.tsv", sep='\t', header=[0, 1, 2], index_col=None)
data.head()
data.columns
s_conc = [3.33, 10, 50, 100]
i_conc = [0, 0.07, 0.15]
plt.figure(figsize=(20, 15))
for i, conc in enumerate(product(i_conc, s_conc)):
plt.subplot(3, 4, i+1)
x = data.loc[:, (f'S = {conc[1]} uM', f'I = {conc[0]} uM', 'time')].dropna()
y = data.loc[:, (f'S = {conc[1]} uM', f'I = {conc[0]} uM', 'O2_consumption')].dropna()
plt.plot(x, y, 'bo')
plt.title(f'S = {conc[1]} uM I = {conc[0]} uM')
plt.xlabel('time')
plt.ylabel('O2 consumption')
plt.tight_layout()
Задание 1. Рассчёт начальной скорости реакции.¶
Воспользуемся методом начальных скоростей и аппроксимируем первые 5 секунд реакции. Для начала посмотрим на эти 5 секунд.
plt.figure(figsize=(20, 15))
for i, conc in enumerate(product(i_conc, s_conc)):
plt.subplot(3, 4, i+1)
time = data.loc[:, (f'S = {conc[1]} uM', f'I = {conc[0]} uM', 'time')].dropna()
cons = data.loc[:, (f'S = {conc[1]} uM', f'I = {conc[0]} uM', 'O2_consumption')].dropna()
plt.plot(time[time < 5], cons[time < 5], 'bo')
plt.title(f'S = {conc[1]} uM I = {conc[0]} uM')
plt.xlabel('time')
plt.ylabel('O2 consumption')
plt.tight_layout()
Данных для каждого эксперимента в первые 5 секунд немного, но они вполне неплохие.
o2_speed = list()
for conc in product(i_conc, s_conc):
time = data.loc[:, (f'S = {conc[1]} uM', f'I = {conc[0]} uM', 'time')].dropna()
cons = data.loc[:, (f'S = {conc[1]} uM', f'I = {conc[0]} uM', 'O2_consumption')].dropna()
slope, intercept = np.polyfit(x=time[time < 5],
y=cons[time < 5],
deg=1)
o2_speed.append([slope, intercept])
o2_speed = np.array(o2_speed, dtype=float)
Для удобства соберём данные в таблицу.
speeds = pd.DataFrame({'speed': o2_speed[:, 0] / 2,
'inhibitor' : [i[0] for i in list(product(i_conc, s_conc))],
'substrate' : [i[1] for i in list(product(i_conc, s_conc))]})
speeds
Посмотрим, как у нас прошли аппроксимации.
plt.figure(figsize=(20, 15))
for i, conc in enumerate(product(i_conc, s_conc)):
plt.subplot(3, 4, i+1)
time = data.loc[:, (f'S = {conc[1]} uM', f'I = {conc[0]} uM', 'time')].dropna()
cons = data.loc[:, (f'S = {conc[1]} uM', f'I = {conc[0]} uM', 'O2_consumption')].dropna()
plt.plot(time[time < 5], cons[time < 5], 'bo')
plt.plot(time[time < 5], time[time < 5] * o2_speed[i, 0] + o2_speed[i, 1], color='red')
plt.title(f'S = {conc[1]} uM I = {conc[0]} uM')
plt.xlabel('time')
plt.ylabel('O2 consumption')
plt.tight_layout()
В целом, вполне хорошо.
Поглядим на зависимость начальной скорости реакции от концентрации субстрата при разных концентрациях ингибитора.
sns.pointplot(data=speeds, x='substrate', y='speed', hue='inhibitor')
Задание 2. Графическое определение типа ингибирование в обратных координатах.¶
speeds['reverse_speed'] = 1 / speeds['speed']
speeds['reverse_substrate'] = 1 / speeds['substrate']
for inh in i_conc:
x = speeds.query(f'inhibitor == {inh}')['reverse_substrate']
y = speeds.query(f'inhibitor == {inh}')['reverse_speed']
plt.plot(x, y, 'o', label=inh)
plt.legend(title='Inhibitor, uM')
plt.xlabel('1 / S')
plt.ylabel('1 / V')
Это всё классно, конечно, но аппроксимировать по 4 точкам — это довольно слабо.
reverse_approx = list()
for inh in i_conc:
x = speeds.query(f'inhibitor == {inh}')['reverse_substrate']
y = speeds.query(f'inhibitor == {inh}')['reverse_speed']
slope, intercept = np.polyfit(x=x,
y=y,
deg=1)
reverse_approx.append((slope, intercept))
for i, inh in enumerate(i_conc):
x = speeds.query(f'inhibitor == {inh}')['reverse_substrate']
y = speeds.query(f'inhibitor == {inh}')['reverse_speed']
plt.plot(x, y, 'o', label=inh)
app_x = np.linspace(-0.3, 0.5, 10)
app_y = app_x * reverse_approx[i][0] + reverse_approx[i][1]
plt.plot(app_x, app_y, color='red')
plt.legend(title='Inhibitor, uM')
plt.xlabel('1 / S')
plt.ylabel('1 / V')
plt.axhline(0, color='black')
plt.axvline(0, color='black')
plt.show()
Судя по графику, прямые не пересекаются в одной точке. Но точно понятно, что прямые не параллельны и что они не пересекаются на оси оординат выше нуля. Вдобавок, точки пересечения близки к оси абсцисс левее нуля, поэтому из анализа графика можно сделать вывод, что здесь наблюдается неконкурентное ингибирование.
Причём, судя по тому, что с ростом концентрации ингибитора прямые вращаются против часовой стрелки, то это действительно ингибирование, а не активация. Иначе говоря, в общей схеме ингибирования, здесь $\alpha = 1$, а $0 \leq \beta \leq 1$.
Задание 3. Определение кинетических параметров реакций.¶
Проведём аппроксимацию начальных скоростей в прямых и двойных обратных координатах для всех указанных в задании механизмов ингибирования и определим точность аппроксимации с помощью MSE.
def competitive_direct(x, vm, km, ki):
s, i = x
return vm * s / (km * (1 + i / ki) + s)
def competitive_reverse(x, vm, km, ki):
s, i = x
return 1 / vm + km * (1 + i / ki) / vm / s
def noncompetitive_direct(x, vm, km, ki):
s, i = x
return vm * s / (1 + i / ki) / (km + s)
def noncompetitive_reverse(x, vm, km, ki):
s, i = x
return (1 + i / ki) / vm + km * (1 + i / ki) / vm / s
def uncompetitive_direct(x, vm, km, ki):
s, i = x
return vm * s / (1 + i / ki) / (km / (1 + i / ki) + s)
def uncompetitive_reverse(x, vm, km, ki):
s, i = x
return (1 + i / ki) / vm + km / vm / s
def mse(x1, x2):
if type(x1) != np.ndarray or type(x2) != np.ndarray:
raise ValueError(f"Arguments are of not supported type: x1 {type(x1)} and x2 {type(x2)}")
return np.square(x1 - x2).mean()
inh_funcs = {func.__name__: func
for func in [competitive_direct,
competitive_reverse,
noncompetitive_direct,
noncompetitive_reverse,
uncompetitive_direct,
uncompetitive_reverse]}
approx_results = list()
approx_keys = ['Vm', 'Km', 'Ki', 'inh_type', 'inh_coord', 'i_conc']
for func_name, func in inh_funcs.items():
inh_type, inh_coord = func_name.split("_")
for inh in i_conc:
approx_data = speeds.query(f'inhibitor == {inh}')
if inh_coord == 'direct':
x = approx_data['substrate']
y = approx_data['speed']
elif inh_coord == 'reverse':
x = approx_data['reverse_substrate']
y = approx_data['reverse_speed']
popt, pcov = curve_fit(func, (x, inh), y, bounds=(0, 'inf'))
approx_results.append(dict(zip(approx_keys, list(popt) + [inh_type, inh_coord, inh])))
inh_results = pd.DataFrame.from_records(approx_results)
inh_results
melted_params = pd.melt(inh_results,
id_vars=['i_conc', 'inh_coord', 'inh_type'],
value_vars=['Ki', 'Km', 'Vm'],
var_name='parameter',
value_name='value')
inh_coord = list(melted_params['inh_coord'].unique())
inh_type = list(melted_params['inh_type'].unique())
print(inh_coord, inh_type)
plt.figure(figsize=(15, 10))
for i, item in enumerate(product(inh_coord, inh_type)):
plt.subplot(2, 3, i + 1)
plot_data = melted_params.query(f'inh_coord == "{item[0]}" and inh_type == "{item[1]}"')
for inh in i_conc:
plot_plot_data = plot_data.query(f"i_conc == {inh}")
plt.plot(plot_plot_data['parameter'], plot_plot_data['value'], 'o', label=inh)
plt.legend(title='Inhibitor, uM')
plt.title(f"{item[1]} inhibition in {item[0]} coordinates")
Мы видим, что в прямых координатах все типы ингибирования согласуются между собой в полученных значениях параметров, в то время как в обратных координатах согласованность пропадает. Если сравнивать прямые и обратные координаты для каждого типа ингибирования, окажется, что:
- У конкурентного ингибирования очень сильно различаются полученные для $K_i$ значения;
- У неконкурентного ингибирования значительно различаются полученные для $K_M$ значения;
- У бесконкурентного ингибирования значительно различаются полученные для $K_M$ и $V_m$ значения.
Но в целом, в пределах 1 порядка данные совпадают только для неконкурентного ингибирования, что ещё раз подтверждает нашу гипотезу о механизме ингибирования.
Давайте посмотрим на $MSE$ скоростей реакций по нашим оценкам кинетических параметров, то есть то, насколько хорошо полученные параметры объясняют реакцию.
mse_results = list()
mse_keys = ['inh_type', 'inh_coord', 'i_conc', 'mse']
for func_name, func in inh_funcs.items():
inh_type, inh_coord = func_name.split("_")
for inh in i_conc:
approx_data = speeds.query(f'inhibitor == {inh}')
if inh_coord == 'direct':
x = approx_data['substrate']
y = approx_data['speed']
elif inh_coord == 'reverse':
x = approx_data['reverse_substrate']
y = approx_data['reverse_speed']
params = inh_results.query(f'inh_coord == "{inh_coord}" and inh_type == "{inh_type}" and i_conc == {inh}')
ki, km, vm = params['Ki'].item(), params['Km'].item(), params['Vm'].item()
y_pred = func((x, inh), vm, km, ki)
error = mse(np.array(y), np.array(y_pred))
mse_results.append(dict(zip(mse_keys, [inh_type, inh_coord, inh, error])))
mse_df = pd.DataFrame.from_records(mse_results)
mse_df
inh_coord = list(melted_params['inh_coord'].unique())
inh_type = list(melted_params['inh_type'].unique())
plt.figure(figsize=(15, 10))
for i, item in enumerate(product(inh_coord, inh_type)):
plt.subplot(2, 3, i + 1)
plot_data = mse_df.query(f'inh_coord == "{item[0]}" and inh_type == "{item[1]}"')
plt.plot(plot_data['i_conc'], plot_data['mse'], 'o')
plt.title(f"{item[1]} inhibition in {item[0]} coordinates")
plt.xlabel("Inhibitor, uM")
plt.ylabel("MSE")
plt.tight_layout()
Удивительно, насколько согласованными оказались значения $MSE$! Я проверил, ошибок в вычислениях нет. Эта одинаковость смущает. Более того, ошибки в прямых координатах маленькие, а в обратных большие.
Давайте исследуем скорости настоящие и предсказанные для разных способов аппроксимации.
speed_results = list()
speed_keys = ['inh_type', 'inh_coord', 'i_conc', 'true_speed', 'predicted_speed']
for func_name, func in inh_funcs.items():
inh_type, inh_coord = func_name.split("_")
for inh in i_conc:
approx_data = speeds.query(f'inhibitor == {inh}')
if inh_coord == 'direct':
x = approx_data['substrate']
y = approx_data['speed']
elif inh_coord == 'reverse':
x = approx_data['reverse_substrate']
y = approx_data['reverse_speed']
params = inh_results.query(f'inh_coord == "{inh_coord}" and inh_type == "{inh_type}" and i_conc == {inh}')
ki, km, vm = params['Ki'].item(), params['Km'].item(), params['Vm'].item()
y_pred = func((x, inh), vm, km, ki)
speed_results.append(dict(zip(speed_keys, [inh_type, inh_coord, inh, y, y_pred])))
speed_compare = pd.DataFrame.from_records(speed_results)
speed_compare
inh_coord = list(melted_params['inh_coord'].unique())
inh_type = list(melted_params['inh_type'].unique())
plt.figure(figsize=(15, 10))
for i, item in enumerate(product(inh_coord, inh_type)):
plt.subplot(2, 3, i + 1)
plot_data = speed_compare.query(f'inh_coord == "{item[0]}" and inh_type == "{item[1]}"')
for inh in i_conc:
plot_plot_data = plot_data.query(f'i_conc == {inh}')
#print(plot_plot_data['true_speed'])
plt.plot(plot_plot_data.loc[:, 'true_speed'].item(), plot_plot_data.loc[:, 'predicted_speed'].item(), 'o', label=inh)
plt.legend(title='Inhibitor, uM')
plt.title(f"{item[1]} inhibition in {item[0]} coordinates")
plt.xlabel("True speed")
plt.ylabel("Predicted speed")
plt.tight_layout()
Что можно понять из этих графиков?
- Скорости предсказываются почти что одинаково между разными типами ингибирования;
- В прямых координатах скорости предсказаны хорошо;
- В обратных координатах скорости предсказаны плохо.
speeds
Как это можно объяснить?
- В обратных координатах концентрации субстрата очень похожи между собой (0.01, 0.02, 0.1 и 0.3) в абсолютных, а скорости отличаются сильнее. По-видимому, из-за этого $MSE$ в обратных координатах большая.
- В прямых координатах концентрации субстрата и скорости хорошо разделимы между собой, поэтому предсказания удаются лучше.
- Ввиду недостатка данных (4 значения субстрата) предсказание параметров (3 шт) не может быть выполнено хорошо. Видимо, поэтому предсказанные скорости и кинетические параметры (в прямых координатах) похожи между разными типами ингибирования.
Теперь определим средние кинетические параметры для разных типов ингибирования.
inh_coords = list(melted_params['inh_coord'].unique())
inh_types = list(melted_params['inh_type'].unique())
mean_records = list()
mean_keys = ['inh_type', 'inh_coord', 'Km', 'Vm', 'Ki']
for inh_coord, inh_type in product(inh_coords, inh_types):
mean_data = inh_results.query(f'inh_coord == "{inh_coord}" and inh_type == "{inh_type}"')
km = mean_data['Km'].mean()
vm = mean_data['Vm'].mean()
ki = mean_data.query("i_conc != 0")['Ki'].mean()
mean_records.append(dict(zip(mean_keys, [inh_type, inh_coord, km, vm, ki])))
mean_df = pd.DataFrame.from_records(mean_records)
mean_df
mean_df.query('inh_coord == "direct" and inh_type == "noncompetitive"')
Задание 4. Аппроксимация в координатах Диксона.¶
Построим зависимость обратной скорости от концентрации ингибитора для каждого значения концентрации субстрата.
for s in s_conc:
plot_data = speeds.query(f"substrate == {s}")
plt.plot(plot_data['inhibitor'], plot_data['reverse_speed'], 'o', label=s)
plt.legend(title="Substrate, uM")
plt.xlabel("Inhibitor, uM")
plt.ylabel("1 / V")
plt.show()
Линейно аппроксимируем зависимости.
dickson_approx = list()
for s in s_conc:
x = speeds.query(f'substrate == {s}')['inhibitor']
y = speeds.query(f'substrate == {s}')['reverse_speed']
slope, intercept = np.polyfit(x=x,
y=y,
deg=1)
dickson_approx.append((slope, intercept))
for i, s in enumerate(s_conc):
x = speeds.query(f'substrate == {s}')['inhibitor']
y = speeds.query(f'substrate == {s}')['reverse_speed']
plt.plot(x, y, 'o', label=s)
app_x = np.linspace(-0.3, 0.15, 10)
app_y = app_x * dickson_approx[i][0] + dickson_approx[i][1]
plt.plot(app_x, app_y, color='red')
plt.legend(title='Substrate, uM')
plt.xlabel('Inhibitor, uM')
plt.ylabel('1 / V')
plt.axhline(0, color='black')
plt.axvline(0, color='black')
plt.show()
Полученные прямые пересекаются на оси абсцисс слева от нуля. Какой это тип ингибирования? Давайте выводить.
В координатах Диксона уравнения для трёх типов ингибирования выглядят следующим образом:
- Конкурентное: $\frac{1}{V} = \frac{1}{V_m} + \frac{K_M}{V_m \cdot S} + \frac{K_M}{V_m \cdot K_i \cdot S} \cdot I$;
- Неконкурентное: $\frac{1}{V} = \frac{1}{V_m} + \frac{K_M}{V_m \cdot S} + (\frac{1}{V_m \cdot K_i} + \frac{K_M}{V_m \cdot K_i \cdot S}) \cdot I$;
- Бесконкурентное: $\frac{1}{V} = \frac{1}{V_m} + \frac{K_M}{V_m \cdot S} + \frac{1}{V_m \cdot K_i} \cdot I$.
Пересечение прямых означает, что $\frac{1}{V}$ не зависит от $S$. Значит, для каждого механизма надо выяснить, при каком значении $I$ $\frac{1}{V}$ не зависит от $S$ и найти это значение $\frac{1}{V}$.
Заметим, что произвольная функция $f(x_1, x_2, ...)$ не зависит от $x_i$ тогда и только тогда, когда $\frac{\partial f}{\partial x_i} = 0$. Таким образом, нам надо взять производные $\frac{\partial}{\partial S}(\frac{1}{V})$ для каждого механизма, приравнять к $0$ и найти искомую точку.
Сделаем это.
- Конкурентное: $\frac{\partial}{\partial S}(\frac{1}{V}) = -\frac{K_M}{V_M \cdot S^2} - \frac{I \cdot K_M}{V_m \cdot K_i \cdot S^2} = 0$, откуда следует, что $I = -K_i$ и $\frac{1}{V} = \frac{1}{V_m}$, то есть точка пересечения во II квадранте;
- Неконкурентное: $\frac{\partial}{\partial S}(\frac{1}{V}) = -\frac{K_M}{V_M \cdot S^2} - \frac{I \cdot K_M}{V_m \cdot K_i \cdot S^2} = 0$, откуда следует, что $I = -K_i$ и $\frac{1}{V} = 0$, то есть точка пересечения на оси абсцисс левее нуля;
- Бесконкурентное: $\frac{\partial}{\partial S}(\frac{1}{V}) = \frac{K_M}{V_m} \neq 0$, то есть точки пересечения нет, а значит, прямые параллельны.
В наших данных мы видим, что у нас неконкурентное ингибирование (согласуется с предыдущими данными).
Теперь определим кинетические параметры в координатах Диксона для трёх механизмов и узнаем точность определения (всё как в задании 3).
def competitive_dickson(x, vm, km, ki):
i, s = x
return 1 / vm + km / vm / s + km * i / vm / ki / s
def noncompetitive_dickson(x, vm, km, ki):
i, s = x
return 1 / vm + km / vm / s + (1 / vm / ki + km / vm / ki / s) * i
def uncompetitive_dickson(x, vm, km, ki):
i, s = x
return 1 / vm + km / vm / s + i / vm / ki
dickson_funcs = {func.__name__: func
for func in [competitive_dickson,
noncompetitive_dickson,
uncompetitive_dickson]}
dickson_records = list()
dickson_keys = ['Vm', 'Km', 'Ki', 'inh_type', 'inh_coord', 's_conc']
for func_name, func in dickson_funcs.items():
inh_type, inh_coord = func_name.split("_")
for s in s_conc:
approx_data = speeds.query(f'substrate == {s}')
x = approx_data['inhibitor']
y = approx_data['reverse_speed']
popt, pcov = curve_fit(func, (x, s), y, bounds=(0, 'inf'))
dickson_records.append(dict(zip(dickson_keys, list(popt) + [inh_type, inh_coord, s])))
dickson_df = pd.DataFrame.from_records(dickson_records)
dickson_df
dickson_df.query('inh_type == "noncompetitive"')
melted_dickson = pd.melt(dickson_df,
id_vars=['s_conc', 'inh_coord', 'inh_type'],
value_vars=['Ki', 'Km', 'Vm'],
var_name='parameter',
value_name='value')
plt.figure(figsize=(15, 5))
for i, item in enumerate(inh_types):
plt.subplot(1, 3, i + 1)
plot_data = melted_dickson.query(f'inh_type == "{item}"')
for s in s_conc:
plot_plot_data = plot_data.query(f"s_conc == {s}")
plt.plot(plot_plot_data['parameter'], plot_plot_data['value'], 'o', label=s)
plt.legend(title='Substrate, uM')
plt.title(f"{item} inhibition in Dickson coordinates")
Наименьший разброс по величинам параметров наблюдается у неконкурентного ингибирования.
Проверим точность предсказаний по полученным параметрам.
mse_records = list()
mse_keys = ['inh_type', 'inh_coord', 's_conc', 'mse']
for func_name, func in dickson_funcs.items():
inh_type, inh_coord = func_name.split("_")
for s in s_conc:
approx_data = speeds.query(f'substrate == {s}')
x = approx_data['inhibitor']
y = approx_data['reverse_speed']
params = dickson_df.query(f'inh_type == "{inh_type}" and s_conc == {s}')
ki, km, vm = params['Ki'].item(), params['Km'].item(), params['Vm'].item()
y_pred = func((x, s), vm, km, ki)
error = mse(np.array(y), np.array(y_pred))
mse_records.append(dict(zip(mse_keys, [inh_type, inh_coord, s, error])))
mse_dickson = pd.DataFrame.from_records(mse_records)
mse_dickson
plt.figure(figsize=(15, 5))
for i, item in enumerate(inh_types):
plt.subplot(1, 3, i + 1)
plot_data = mse_dickson.query(f'inh_type == "{item}"')
plt.plot(plot_data['s_conc'], plot_data['mse'], 'o')
plt.title(f"{item} inhibition in Dickson coordinates")
plt.xlabel("Substrate, uM")
plt.ylabel("MSE")
Для всех механизмов есть небольшая ошибка при малой концентрации субстрата, а так она очень мала. По ошибке нельзя ничего сказать о том, какой механизм лучше.
Посчитаем средние значения параметров.
mean_records = list()
mean_keys = ['inh_type', 'inh_coord', 'Km', 'Vm', 'Ki']
for inh_type in inh_types:
mean_data = dickson_df.query(f'inh_type == "{inh_type}"')
km = mean_data['Km'].mean()
vm = mean_data['Vm'].mean()
ki = mean_data['Ki'].mean()
mean_records.append(dict(zip(mean_keys, [inh_type, inh_coord, km, vm, ki])))
mean_dickson = pd.DataFrame.from_records(mean_records)
mean_dickson
mean_dickson.query('inh_type == "noncompetitive"')
Сравним полученные средние значения параметров с тем, что получили в прямых и двойных обратных координатах.
all_means = pd.concat([mean_dickson, mean_df], axis=0)
all_means
melted_means = pd.melt(all_means, id_vars = ['inh_coord', 'inh_type'], value_vars=['Ki', 'Km', 'Vm'], var_name='parameter')
plt.figure(figsize=(15, 5))
for i, inh_coord in enumerate(['dickson', 'direct', 'reverse']):
plt.subplot(1, 3, i + 1)
sns.scatterplot(data=melted_means.query(f'inh_coord == "{inh_coord}"'), x='parameter', y='value', hue='inh_type')
plt.ylim(0, 8)
plt.title(inh_coord)
plt.tight_layout()
Кинетические параметры более-менее согласуются для трёх механизмов в координатах Диксона и прямых координатах. При этом, наибольшие проблемы вызвала $K_M$ — то у неё большой разброс, то экстремально маленькие значения (в обратных координатах).
Из полученных результатов нельзя сказать, какие координаты лучше использовать — Диксона или прямые. В жизни это зависит от задачи: если у нас много разных концентраций субстрата, то логично использовать прямые координаты, а если много разных концентраций ингибитора, то Диксона.
Задание 5. Совместная аппроксимация.¶
Будем аппроксимировать скорости реакций в прямых координатах и координатах Диксона сразу по всем концентрациям субстрата и ингибитора.
joint_funcs = {func.__name__: func
for func in [competitive_direct,
noncompetitive_direct,
uncompetitive_direct,
competitive_dickson,
noncompetitive_dickson,
uncompetitive_dickson]}
joint_records = list()
joint_keys = ['Vm', 'Km', 'Ki', 'inh_type', 'inh_coord']
for func_name, func in joint_funcs.items():
inh_type, inh_coord = func_name.split("_")
x_s = speeds['substrate']
i_s = speeds['inhibitor']
if inh_coord == 'direct':
y = speeds['speed']
x = [x_s, i_s]
elif inh_coord == 'dickson':
y = speeds['reverse_speed']
x = [i_s, x_s]
popt, pcov = curve_fit(func, x, y, bounds=(0, 'inf'))
joint_records.append(dict(zip(joint_keys, list(popt) + [inh_type, inh_coord])))
joint_df = pd.DataFrame.from_records(joint_records)
joint_df
melted_joint = pd.melt(joint_df, id_vars=['inh_coord', 'inh_type'], value_vars=['Ki', 'Km', 'Vm'], var_name='parameter')
Посмотрим на то, как разные механизмы дают параметры в одних и тех же координатах.
plt.figure(figsize=(10, 5))
for i, item in enumerate(['dickson', 'direct']):
plt.subplot(1, 2, i + 1)
plot_data = melted_joint.query(f'inh_coord == "{item}"')
sns.scatterplot(data=plot_data, x='parameter', y='value', hue='inh_type')
plt.title(f"Inhibition in {item} coordinates")
Как и в других наблюдениях, $K_i$ и $V_m$ похожи, а $K_M$ различаются, притом весьма сильно, как между координатами, так и внутри них.
Посмотрим на то, как один и тот же механизм работает в разных координатах.
plt.figure(figsize=(15, 5))
for i, item in enumerate(inh_types):
plt.subplot(1, 3, i + 1)
plot_data = melted_joint.query(f'inh_type == "{item}"')
sns.scatterplot(data=plot_data, x='parameter', y='value', hue='inh_coord')
plt.title(f"{item} inhibiton")
Как и в других наблюдениях, $K_i$ и $V_m$ похожи, а $K_M$ различаются для одного и того же механизма в разных координатах. При этом, для неконкурентного ингибирования различие в $K_i$ и $V_m$ минимально.
print("MSE\ttype")
for item in inh_types:
compare_data = joint_df.query(f'inh_type == "{item}"')
direct_data = compare_data.query(f'inh_coord == "direct"').loc[:, ['Ki', 'Km', 'Vm']].T
dickson_data = compare_data.query(f'inh_coord == "dickson"').loc[:, ['Ki', 'Km', 'Vm']].T
error = mse(np.array(direct_data), np.array(dickson_data))
print(f"{error:.4f}\t{item}")
Оказывается, что у конкурентного ингибирования $MSE$ различий между координатами меньше, чем у неконкурентного. Это достигнуто меньшей разницей в значениях $K_M$. И вновь вывод в том, что величина $MSE$ не может быть критерием отбора механизмов.
plt.figure(figsize=(15, 10))
for i, item in enumerate(inh_types):
plt.subplot(2, 3, i + 1)
plot_data = melted_joint.query(f'inh_type == "{item}"')
sns.scatterplot(data=plot_data, x='parameter', y='value', hue='inh_coord')
plt.title(f"{item} inhibiton by joint approximation")
for i, item in enumerate(inh_types):
plt.subplot(2, 3, i + 4)
plot_data = melted_means.query(f'inh_type == "{item}"')
sns.scatterplot(data=plot_data, x='parameter', y='value', hue='inh_coord')
plt.ylim(0, 20)
plt.title(f"{item} inhibiton by mean data")
Параметры, полученные для одного и того же механизма разными способами, сходятся между собой в разных координатах.
Давайте порисуем в 3D.
from mpl_toolkits.mplot3d import Axes3D
Нарисуем проволочную поверхность того, как предсказывают поведение системы разные типы ингибирования в прямых координатах и координатах Диксона.
fig = plt.figure(figsize=(15, 10))
for i, (func_name, func) in enumerate(joint_funcs.items()):
inh_type, inh_coord = func_name.split("_")
substrate = s_conc
inhibitor = i_conc
xs, ys = np.meshgrid(substrate, inhibitor)
ki, km, vm, _, _ = joint_df.query(f'inh_coord == "{inh_coord}" and inh_type == "{inh_type}"').values.squeeze()
if inh_coord == 'direct':
zs = func((xs, ys), vm, km, ki)
elif inh_coord == 'dickson':
zs = func((ys, xs), vm, km, ki)
ax = fig.add_subplot(2, 3, i+1, projection='3d')
ax.plot_wireframe(xs, ys, zs)
ax.scatter(speeds['substrate'],
speeds['inhibitor'],
speeds['speed'] if inh_coord == 'direct' else speeds['reverse_speed'],
color='red')
plt.title(f"{inh_coord} {inh_type}")
ax.set_xlabel("substrate, uM")
ax.set_ylabel("inhibitor, uM")
ax.set_zlabel("V" if inh_coord == 'direct' else "1 / V")
plt.tight_layout()
Неконкурентное ингибирование наилучшим образом описывает поведение системы в данных координатах.
Посмотрим на предсказываемое состояние системы в других точках.
fig = plt.figure(figsize=(15, 10))
for i, (func_name, func) in enumerate(joint_funcs.items()):
inh_type, inh_coord = func_name.split("_")
substrate = np.logspace(0, 2, 20)
inhibitor = np.linspace(0, 0.2, 20)
xs, ys = np.meshgrid(substrate, inhibitor)
ki, km, vm, _, _ = joint_df.query(f'inh_coord == "{inh_coord}" and inh_type == "{inh_type}"').values.squeeze()
if inh_coord == 'direct':
zs = func((xs, ys), vm, km, ki)
elif inh_coord == 'dickson':
zs = func((ys, xs), vm, km, ki)
ax = fig.add_subplot(2, 3, i+1, projection='3d')
ax.plot_surface(xs, ys, zs, color='#1F77B4DD')
ax.scatter(speeds['substrate'],
speeds['inhibitor'],
speeds['speed'] if inh_coord == 'direct' else speeds['reverse_speed'],
color='red')
plt.title(f"{inh_coord} {inh_type}")
ax.set_xlabel("substrate, uM")
ax.set_ylabel("inhibitor, uM")
ax.set_zlabel("V" if inh_coord == 'direct' else "1 / V")
plt.tight_layout()
Зачем? Просто красиво. Можно увидеть, что неконкурентное и бесконкурентное ингибирования очень аккуратно описывают систему в прямых координатах (но неконкуретное аккуратнее). В координатах Диксона лучше себя ведёт неконкурентное, правильнее описывая всю систему в целом и поведение около нулевых концентраций в частности.
Задание 6. Рассчёт общего случая.¶
Рассчитаем параметры системы в общем случае.
def general_direct(x, vm, km, ki, a, b):
s, i = x
return vm * s * (1 + b * i / a / ki) / (km * (1 + i / ki) + s * (1 + i / a / ki))
popt, pcov = curve_fit(general_direct,
[speeds['substrate'],
speeds['inhibitor']],
speeds['speed'],
bounds=(0, 'inf'))
vm, km, ki, a, b = list(popt.squeeze())
print("Vm\tKm\tKi\ta\tb")
print(f"{vm:.4f}\t{km:.4f}\t{ki:.4f}\t{a:.4f}\t{b}")
fig = plt.figure(figsize=(5, 5))
substrate = s_conc
inhibitor = i_conc
xs, ys = np.meshgrid(substrate, inhibitor)
zs = general_direct((xs, ys), vm, km, ki, a, b)
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(xs, ys, zs, color='#1F77B4DD')
ax.scatter(speeds['substrate'],
speeds['inhibitor'],
speeds['speed'],
color='red')
plt.title("General inhibiton")
ax.set_xlabel("Substrate, uM")
ax.set_ylabel("Inhibitor, uM")
ax.set_zlabel("V, uM/s")
plt.tight_layout()
Нарисовало практически идеально.
fig = plt.figure(figsize=(5, 5))
substrate = np.logspace(0, 2, 20)
inhibitor = np.linspace(0, 0.2, 20)
xs, ys = np.meshgrid(substrate, inhibitor)
zs = general_direct((xs, ys), vm, km, ki, a, b)
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(xs, ys, zs, color='#1F77B4DD')
ax.scatter(speeds['substrate'],
speeds['inhibitor'],
speeds['speed'],
color='red')
plt.title("General inhibiton")
ax.set_xlabel("Substrate, uM")
ax.set_ylabel("Inhibitor, uM")
ax.set_zlabel("V, uM/s")
plt.tight_layout()
И весьма красиво :)
Посмотрим, какие скорости выдаст $K_M$ как концентрация субстрата.
half_v = general_direct((km, np.array(i_conc)), vm, km, ki, a, b)
half_v
sns.pointplot(data=speeds, x='substrate', y='speed', hue='inhibitor')
plt.axvline(km / 10, color='red') # не знаю, почему такая проблема.
plt.axhline(half_v[0], color='navy')
plt.axhline(half_v[1], color='orange')
plt.axhline(half_v[2], color='green')
plt.ylim(0, 1)
plt.show()
Всё сходится со здравым смыслом.