Python、scipyのcurve_fitを使ってミカエリスメンテン式へのfittingを行う

Pythonのscipyを使ってcurve_fitを利用したミカエリスメンテン式へのfittingを行う。
Pandasを利用してデータを読み込み、seabornを用いてグラフを出力します。
この方法を用いてGoogle Colaboratoryでも描画できます。

https://ja.wikipedia.org/wiki/%E3%83%9F%E3%82%AB%E3%82%A8%E3%83%AA%E3%82%B9%E3%83%BB%E3%83%A1%E3%83%B3%E3%83%86%E3%83%B3%E5%BC%8F



<h1>モジュールのインポート</h1>
まず、必要なモジュールをimportする。

```
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit
import seaborn as sns
```

<h1>データの準備</h1>
Pandas形式のデータを作成する。
直接コードに書き込む場合は以下の方法を用いると良い。
1, 2, 3, 4 種類の濃度のデータを3回解析したデータを用いるとする。

```
#濃度のリスト
x = [1,2,3,4,
     1,2,3,4,
     1,2,3,4,]
#活性のリスト
y = [0.0010,0.0020,0.0025,0.0030,
     0.0015,0.0025,0.0027,0.0032,
     0.0008,0.0017,0.0023,0.0028]
data_1 = pd.DataFrame({'concentration': x, 'activity': y,})
```
csvファイルの読み込みでもよい。
以下のような中身のcsvファイルに任意の名前をつけてcsvファイルを読み込むと良い。

```
concentration,activity
1,0.0010
2,0.0020
3,0.0025
4,0.0030
1,0.0015
2,0.0025
3,0.0027
4,0.0032
1,0.0008
2,0.0017
3,0.0023
4,0.0028
```

以下のコマンドでcsvファイルを読み込む。
google colaboratoryを利用する場合は

```
#Google colaboratoryにファイルをアップロードする場合
uploaded = files.upload()

#csvファイルの読み込み
data = pd.read_csv("filename.csv")
```

printを用いてデータリストが正しく作成されていることを確認する。
濃度と活性のデータが1対１で対応していることを確認します。
一つの濃度につき、3種のデータがあります。

```
print(data_1)
#以下出力
    concentration  activity
0               1    0.0010
1               2    0.0020
2               3    0.0025
3               4    0.0030
4               1    0.0015
5               2    0.0025
6               3    0.0027
7               4    0.0032
8               1    0.0008
9               2    0.0017
10              3    0.0023
11              4    0.0028
```

試しにこれらのデータを散布図として出力し、傾向を確認する。
```
sns.scatterplot(data = data_1, x = "concentration", y = "activity")
```
以下のようなグラフが出力されます。

![ダウンロード.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1825372/6ede9dd3-1202-c8cb-64db-62ef8844e6c5.png)



<h1>双曲線にfitting</h1>
まず、fittingに用いる関数を用意します。

```
#ミカエリスメンテン式
#k1がVmax, k2がKmを表しています。
def MM_Eq(x,k1,k2): 
    return  k1*x/(x+k2)
```

次にcurve fitを用いてfittingを行います。

```
param, cov = curve_fit(MM_Eq,data_1["concentration"],data_1["activity"])
print(param)
print(cov)
```
paramにfitting後のk1とk2が格納されます。
k1がVmax, k2がKmに対応します。
covには共分散が出力されます。

次に出力された近似曲線を描画するためのdefを作ります。
この方法により、fittingした結果得られた関数の理想値のリストを作成します。
ここでは100個の点を作成し、曲線の描画に用います。
リスト作成のためのdefをまず用意します。
```
#近似式に該当した曲線の出力
#曲線を描画するために用いる点の数を指定します。
RES = 100

#maxには濃度の最大値を入力します。
#k1, k2は先ほどの値と同様です。
def graph_plot(max, k1, k2): 
  limit = max + max/10
  x = []
  y = []
  lst = range(0,RES + 1)
  for j in lst:
    i = max/RES * j
    x.append(i)
    y.append(MM_Eq(i,k1,k2))
  return x, y
```

作成したdefを用いて、実際に点のデータを作成します。
さらにpandas形式に変換し格納します。
data["concentration"].max()は基質濃度の最大値を返します。
param[0],param[1]はそれぞれ出力されたk1とk2を表します。
```
xs, ys = graph_plot(data_1["concentration"].max(),param[0],param[1])
curve = pd.DataFrame({'concentration': xs, 'activity': ys,})
```

pandasのリストが正しくできていることを確認します
最大濃度が4、点の数を100 (RES)と指定したため、0.04 ごとにfittingした関数の結果が出力されたリストが得られます。
```
print(curve)
#以下が出力
     concentration  activity
0             0.00  0.000000
1             0.04  0.000057
2             0.08  0.000113
3             0.12  0.000167
4             0.16  0.000221
..             ...       ...
96            3.84  0.002932
97            3.88  0.002948
98            3.92  0.002964
99            3.96  0.002980
100           4.00  0.002995
```

試しに散布図を描画します。

```
sns.scatterplot(data = curve, x = "concentration", y = "activity", )
```
![ダウンロード (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1825372/36c58c79-0752-6bc4-d922-a91a3be79e85.png)
散布図なので、点が並んでいますが、最終的には折れ線グラフを用いた方がいいと思われます。

<h1>グラフの描画</h1>
元、データとfittingした関数を重ね合わせた図を出力します。

```
#入力データの描画
sns.scatterplot(data = data_1, x = "concentration", y = "activity")
#fittingした関数の描画
sns.lineplot(data = curve, x = "concentration", y = "activity")
#pdf形式でファイルを保存する
plt.savefig('filename.pdf')
```

以下のような画像がダウンロードできます。
![ダウンロード (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1825372/cbeae3bf-b5a0-c6ef-99fd-d59450abd4be.png)


平均と標準誤差をプロットしたい場合

```
#濃度の抽出
data_conc = data_1.concentration.drop_duplicates()

#活性の平均の計算
data_ave = data_1.groupby('concentration').mean()

#活性の平均の標準誤差の計算
data_se = data_1.groupby('concentration').sem()

#図形の描画
plt.errorbar(data_conc.values,data_ave.values, 
             yerr = data_se.values, capsize=5, 
             fmt='o', markersize=10, ecolor='black', 
             markeredgecolor = "black", color='w')
plt.plot(curve['concentration'],curve['activity'], color = 'black')
plt.savefig('filename.pdf')

```
この場合、以下のような図が出力されます。

![ダウンロード-10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1825372/321f43ee-e02f-3570-51f7-e3d938da2c98.png)