第1週、1日目：モデルの種類

Neuromatch Academyによる提供

コンテンツ制作者： マット・ラポート、バイロン・ギャルブレス、コンラッド・コーディング

コンテンツレビュワー： ダリン・グオ、アイシュワリヤ・バルワニ、マディネ・サルヴェスタニ、マリヤム・ヴァジリ・パシカム、マイケル・ワスコム、エラ・バティ

制作編集者： ガガナ・ビー、スピロス・チャブリス

ここで使用されているデータの一部は、Steinmetz et al. (2019)によって共有されたものであることをご承知おきください。

チュートリアルの目標 ​

このチュートリアルは、神経データを理解するために使用される異なる種類のモデルに関する3部シリーズの1つ目です。このチュートリアルでは、「What」モデルを探求します。これは、データを記述するために使用されます。データの特徴を理解するために、異なる方法で視覚化します。そして、それを簡単な数学モデルと比較します。具体的には、次のことを行います：

数百のニューロンからのスパイク活動を含むデータセットを読み込み、その構造を理解します。 ポップレーション全体のスパイク活動の特性を視覚化するためのプロットを作成します。 単一のニューロンの「インタースパイクインターバル」（ISI）の分布を計算します。 この分布の形状に関するいくつかの形式的モデルを考慮し、手動でデータに適合させます。

TODO: コードを転記

Copy code スライドをダウンロードしたい場合は、こちらから: https://osf.io/download/6dxwe/

TODO: スライドを転記

セットアップ ​

フィードバックガジェットのインストールとインポート ​

TODO: ここにコードを転記

Pythonでは、関数を使用する前にライブラリを明示的に「import」する必要があります。私たちは常に各ノートブックまたはスクリプトの冒頭にインポートを指定します。

javascript Copy code import numpy as np import matplotlib.pyplot as plt

チュートリアルのノートブックは、通常、デフォルトでは表示されないいくつかのセットアップ手順で始まります。

重要： コードが非表示になっていても、ノートブックの残りの部分が正常に動作するように、コードを実行する必要があります。各セルを進めるには、左上隅にある再生ボタンを押すか、キーボードショートカット（MacではCmd-Return、その他の場合はCtrl-Enter）を使用します。番号が括弧内に表示され（例：[3]）、そのセルが実行されたことと、それが実行された順序が示されます。

各セルの中身を確認したい場合は、ダブルクリックして展開することができます。展開されたら、エディタの右側の白いスペースをダブルクリックして再び折りたたむことができます。

“What”モデル ​

ビデオ1: ​

ビデオは以下のリンクからご覧いただけます：https://youtube.com/watch?v=KgqR_jbjMQg

セクション1：Steinmetzデータセットの探索 ​

このチュートリアルでは、神経科学のデータセットの構造を探索します。

Steinmetz et al. (2019)の研究の一部のデータを考慮します。この研究では、Neuropixelsプローブがマウスの脳にインプラントされました。各プローブの長さに沿って数百の電極で電位が測定されました。各電極の測定値は、近くのスパイクするニューロンによる電場の局所的な変動を捉えています。スパイクのソートアルゴリズムが使用され、スパイクの時刻が推定され、共通の起源に基づいてスパイクがクラスタリングされました：ソートされたスパイクの単一のクラスタは、単一のニューロンに因果関係を持たせています。

特に、スパイクの時刻とニューロンの割り当てがロードされ、前述のセットアップでspike_timesに割り当てられました。

通常、データセットにはその構造に関する情報が付属しています。ただし、この情報が不完全な場合もあります。また、データの興味を持つ部分を作成するためにいくつかの変換や「前処理」を適用することがありますが、これは状況によっては部分的に文書化されていない場合があります。いずれにしても、利用可能なツールを使用してデータ構造の不慣れな側面を調査できるようにすることが重要です。

それでは、データがどのように見えるかを見てみましょう

セクション1.1: spike_timesを使ってウォームアップ ​

変数のPythonの型は何ですか？

type(spike_times)

numpy.ndarray

numpy.ndarrayが表示されるはずです。これは通常のNumPy配列を意味します。

エラーメッセージが表示される場合は、おそらくノートブックの上部にあるセットアップセルを実行していないことを意味します。その場合は、セットアップセルを実行してください。

すべてが正常に動作している場合、次の質問をデータセットについて行うことができます：その形は何ですか？

spike_times.shape

(734,)

1次元の要素が734個あり、他の次元はありません。最初のエントリのPythonの型は何であり、その形は何ですか？

idx = 0
print(
  type(spike_times[idx]),
  spike_times[idx].shape,
  sep="\n",
)

<class 'numpy.ndarray'>
(826,)

これも1次元のNumPy配列です！なぜこれがspike_timesの形状に2番目の次元として表示されなかったのでしょうか？つまり、なぜspike_times.shape == (734, 826)ではないのでしょうか？

調査するために、別のエントリをチェックしてみましょう。

idx = 321
print(
  type(spike_times[idx]),
  spike_times[idx].shape,
  sep="\n",
)

<class 'numpy.ndarray'>
(9723,)

これも1次元のNumPy配列であり、形状が異なります。これらの配列内の値のNumPy型と、最初のいくつかの要素をチェックすると、浮動小数点数（別のnp.ndarrayのレベルではなく）から構成されていることがわかります。

i_neurons = [0, 321]
i_print = slice(0, 5)

for i in i_neurons:
  print(
    "ニューロン {}:".format(i),
    spike_times[i].dtype,
    spike_times[i][i_print],
    "\n",
    sep="\n"
  )

ニューロン 0:
float32
[ 0.8149   14.822467 24.9646   25.1436   38.8709  ]

ニューロン 321:
float32
[1.0698667 1.1536334 1.2403667 1.7072    1.799    ]

今回はPythonの変数の型ではなく、NumPyのdtypeをチェックしています。これらの2つの配列は、32ビットの精度を持つ浮動小数点数（「floats」）を含んでいます。

基本的な構造は次のようになります：

• spike_timesは1次元であり、そのエントリはNumPy配列です。その長さはニューロンの数（734）です。インデックスを使用して、ニューロンの一部を選択します。
• spike_times内の配列も1次元であり、単一のニューロンに対応します。そのエントリは浮動小数点数であり、その長さはそのニューロンに帰属するスパイクの数です。インデックスを使用して、そのニューロンのスパイクタイムの一部を選択します。

視覚的には、データ構造は次のように見えると考えることができます：

| . . . . . |
| . . . . . . . . |
| . . . |
| . . . . . . . |

次に進む前に、データセット内のニューロンの数とニューロンごとのスパイク数を計算して保存します。

n_neurons = len(spike_times)
total_spikes_per_neuron = [len(spike_times_i) for spike_times_i in spike_times]

print(f"Number of neurons: {n_neurons}")
print(f"Number of spikes for first five neurons: {total_spikes_per_neuron[:5]}")

Number of neurons: 734
Number of spikes for first five neurons: [826, 2818, 3953, 646, 1115]

ビデオ2: ​

ビデオは以下のリンクからご覧いただけます：https://youtube.com/watch?v=oHwYWUI_o1U

セクション1.2：暖まってきました：合計スパイク数のカウントとプロット ​

これまでにも見てきたように、録音全体のスパイク数はニューロンごとに異なります。より一般的には、特定の期間において、一部のニューロンは他のニューロンよりもより多くスパイクする傾向があります。データセット内のすべてのニューロンにおけるスパイクの分布を探索してみましょう。

ほとんどのニューロンは平均に比べて「活発」なのか、「静か」なのか？これを確認するために、合計スパイク数を基準に一定の幅のビンを定義し、各ビンに含まれるニューロンの数をカウントします。これは「ヒストグラム」として知られています。

matplotlib関数plt.histを使用してヒストグラムをプロットすることができます。計算するだけであれば、numpy関数np.histogramを代わりに使用することができます。

plt.hist(total_spikes_per_neuron, bins=50, histtype="stepfilled") plt.xlabel("Total spikes per neuron") plt.ylabel("Number of neurons");

平均スパイク数よりも少ないスパイク数を持つニューロンの割合を見てみましょう：

mean_spike_count = np.mean(total_spikes_per_neuron)
frac_below_mean = (total_spikes_per_neuron < mean_spike_count).mean()
print(f"ニューロンの{frac_below_mean:2.1%}が平均以下です")

ニューロンの68.0%が平均以下です

また、ヒストグラムプロットに平均スパイク数を追加して表示できます：

plt.hist(total_spikes_per_neuron, bins=50, histtype="stepfilled")
plt.xlabel("ニューロンごとの合計スパイク数")
plt.ylabel("ニューロンの数")
plt.axvline(mean_spike_count, color="orange", label="平均のニューロン")
plt.legend();

このグラフから、多くのニューロンが平均に比べて比較的「静か」であり、一部のニューロンが非常に「活発」であることがわかります。これらのニューロンはより多くスパイクして大きな数を達成する必要があります。

コーディング演習1.2: 平均と中央値のニューロンの比較 ​

もし平均のニューロンが人口の68%よりも活動的である場合、平均のニューロンと中央値のニューロンの関係について何が言えるでしょうか？

演習の目的: 上記のプロットを再現し、中央値のニューロンを追加します。

#################################################################################

TODO for students: ​

Fill out function and remove ​

raise NotImplementedError("Student exercise: complete histogram plotting with median") #################################################################################

Compute median spike count ​

median_spike_count = ... # Hint: Try the function np.median

Visualize median, mean, and histogram ​

plt.hist(..., bins=50, histtype="stepfilled") plt.axvline(..., color="limegreen", label="Median neuron") plt.axvline(mean_spike_count, color="orange", label="Mean neuron") plt.xlabel("Total spikes per neuron") plt.ylabel("Number of neurons") plt.legend()

Click for solution

ボーナス: 中央値は50パーセンタイルです。他のパーセンタイルはどうでしょうか？ヒストグラムに四分位範囲を表示できますか？ ​

フィードバックの提出 ​

セクション2：ニューロンのスパイク活動の視覚化 ​

チュートリアルの開始からここまでの所要時間の目安：15分

セクション2.1：データのサブセットの取得 ​

これからスパイクのトレインを視覚化します。記録が長いため、まず短い時間区間を定義し、この区間内のスパイクのみを視覚化します。これを行うためのヘルパー関数restrict_spike_timesを定義しました。この関数にhelp()を呼び出すと、関数について少し説明が表示されます。

ヘルパー関数restrict_spike_timesのセルを実行してください。

help(restrict_spike_times)

Help on function restrict_spike_times in module __main__:

restrict_spike_times(spike_times, interval)
    Given a spike_time dataset, restrict to spikes within given interval.

    Args:
      spike_times (sequence of np.ndarray): List or array of arrays,
        each inner array has spike times for a single neuron.
      interval (tuple): Min, max time values; keep min <= t < max.

    Returns:
      np.ndarray: like `spike_times`, but only within `interval`

t_interval = (5, 15)  # units are seconds after start of recording
interval_spike_times = restrict_spike_times(spike_times, t_interval)

これは代表的な区間でしょうか？この区間に含まれる合計スパイク数の割合はどれくらいですか？

original_counts = sum([len(spikes) for spikes in spike_times])
interval_counts = sum([len(spikes) for spikes in interval_spike_times])
frac_interval_spikes = interval_counts / original_counts
print(f"合計スパイク数の{frac_interval_spikes:.2%}がこの区間に含まれています")

合計スパイク数の0.33%がこの区間に含まれています

この値は、区間の期間と実験の期間との比率と比較してどうでしょうか？（この区間の合計時間は、全体の時間に対してどれくらいの割合ですか？）

実験の期間を近似的に計算するために、データセット全体での最小スパイク時刻と最大スパイク時刻を取得します。これを行うために、すべてのニューロンを1つの配列に「連結」し、np.ptp（ピークからピーク）を使用して最大値と最小値の差を取得します：

spike_times_flat = np.concatenate(spike_times)
experiment_duration = np.ptp(spike_times_flat)
interval_duration = t_interval[1] - t_interval[0]

frac_interval_time = interval_duration / experiment_duration
print(f"合計時間の{frac_interval_time:.2%}がこの区間に含まれています")

合計時間の0.37%がこの区間に含まれています

これらの2つの値-合計スパイク数の割合と合計時間の割合-は類似しています。これは、この区間におけるニューロンの平均スパイクレートが、全体の記録と比較してあまり異ならないことを示唆しています。

セクション2.2：スパイクトレインとラスターのプロット ​

代表的なサブセットを取得したので、matplotlibのplt.eventplot関数を使ってスパイクをプロットしましょう。まずは単一のニューロンを見てみます：

neuron_idx = 1
plt.eventplot(interval_spike_times[neuron_idx], color=".2")
plt.xlabel("時間（秒）")
plt.yticks([]);

複数のニューロンをプロットすることもできます。以下は3つのニューロンのプロット例です：

neuron_idx = [1, 11, 51]
plt.eventplot(interval_spike_times[neuron_idx], color=".2")
plt.xlabel("時間（秒）")
plt.yticks([]);

これは「ラスタープロット」と呼ばれるもので、各ニューロンのスパイクが異なる行に表示されます。

多数のニューロンをプロットすると、集団内の特性が分かります。記録されたすべてのニューロンのうち、5番目のニューロンごとにプロットしてみましょう：

neuron_idx = np.arange(0, len(spike_times), 5)
plt.eventplot(interval_spike_times[neuron_idx], color=".2")
plt.xlabel("時間（秒）")
plt.yticks([]);

質問: このプロットの情報は、上記で見た合計スパイク数のヒストグラムとどのように関連していますか？

Click for solution

セクション3: インタースパイク間隔とその分布 ​

チュートリアル開始からこの地点までの推定所要時間：25分

spike_timesに含まれる各ニューロンのスパイク時刻の順序付き配列がわかりました（先ほど可視化しました）。次に何を尋ねることができるでしょうか？

科学的な質問は既存のモデルによって形成されます。したがって、このデータに関する質問を導くためにどのような知識を持っていますか？

ニューロンの発火には物理的な制約があることを知っています。発火にはエネルギーがかかり、そのエネルギーはニューロンの細胞内機構が有限の速度で取得できるものです。したがって、ニューロンには不応期があります：ニューロンはその代謝プロセスがサポートできる限りの速さでしか発火できず、同じニューロンの連続したスパイクの間には最小の遅延があります。

より一般的には、「ニューロンは次にどれだけの時間を待つのか？」や「ニューロンが最も長く待つ時間はどれだけか？」といった質問により直接的に対処するために、スパイク時刻を他の何かに変換できるかもしれません。

ここでは、インタースパイク間隔（interspike intervals：ISIs）を考えることができます。これは、同じニューロンの連続するスパイクの間の時間差です。

Exercise 3：単一のニューロンのISIの分布をプロットする ​

エクササイズの目的: スパイク数のヒストグラムと同様に、データセット内のニューロンの1つに対するISIの分布を示すヒストグラムを作成します。

これを3つのステップで行います：

1. 1つのニューロンのスパイク時刻を抽出します
2. ISIs（スパイク間の時間量、または同じスパイク時刻の隣接する2つのスパイク間の差）を計算します
3. 個々のISIの配列でヒストグラムをプロットします

def compute_single_neuron_isis(spike_times, neuron_idx):
  """Compute a vector of ISIs for a single neuron given spike times.

  Args:
    spike_times (list of 1D arrays): Spike time dataset, with the first
      dimension corresponding to different neurons.
    neuron_idx (int): Index of the unit to compute ISIs for.

  Returns:
    isis (1D array): Duration of time between each spike from one neuron.
  """
  #############################################################################
  # Students: Fill in missing code (...) and comment or remove the next line
  raise NotImplementedError("Exercise: compute single neuron ISIs")
  #############################################################################

  # Extract the spike times for the specified neuron
  single_neuron_spikes = ...

  # Compute the ISIs for this set of spikes
  # Hint: the function np.diff computes discrete differences along an array
  isis = ...

  return isis

# Compute ISIs
single_neuron_isis = compute_single_neuron_isis(spike_times, neuron_idx=283)

# Visualize ISIs
plot_isis(single_neuron_isis)

Click for solution

Example output:

一般的に、ISIの短い間隔が優勢であり、ISIが増加するにつれてカウントが急速に減少しています（比較的滑らかに）。ただし、分布の最大値（8〜11 ms）以下のISIでは、カウントも急速に

ゼロに減少しています。これらの非常に低いISIが存在しないことは、不応期仮説と一致しています：ニューロンはこのISIの領域を埋めるほど速く発火することはできません。

他のニューロンの分布も確認してみてください。分布のさまざまな特徴を解決するには、plt.histの呼び出しでbinsの値を調整する必要があるかもしれません。少なすぎるビンを使用すると興味深い詳細がスムーズになるかもしれませんが、ビンが多すぎるとランダムな変動が支配するようになるかもしれません。

また、分布の形状をより詳細に観察するために、plt.histのrange引数を使用して範囲を制限することもできます。 ヒント: plt.histはrange引数を受け取ります。

セクション4: インタースパイク間隔（ISI）の分布の関数形は何ですか？ ​

チュートリアル開始からこの地点までの推定所要時間：35分

Video 4: ISIの分布 ​

• Youtube
• Bilibili

Video available at https://youtube.com/watch?v=DHhM80MOTe8

フィードバックを提出する ​

ISIのヒストグラムは、最大値以上では連続的で単調減少する関数に従うようです。この関数は明らかに非線形です。これは単一の関数族に属する可能性がありますか？

生理現象を説明するために数学的な関数を使用するアイデアを示すために、以下のような異なる関数形を定義してみましょう：指数関数、逆数、および線形関数。

def exponential(xs, scale, rate, x0):
  """A simple parameterized exponential function, applied element-wise.

  Args:
    xs (np.ndarray or float): Input(s) to the function.
    scale (float): Linear scaling factor.
    rate (float): Exponential growth (positive) or decay (negative) rate.
    x0 (float): Horizontal offset.

  """
  ys = scale * np.exp(rate * (xs - x0))
  return ys

def inverse(xs, scale, x0):
  """A simple parameterized inverse function (`1/x`), applied element-wise.

  Args:
    xs (np.ndarray or float): Input(s) to the function.
    scale (float): Linear scaling factor.
    x0 (float): Horizontal offset.

  """
  ys = scale / (xs - x0)
  return ys

def linear(xs, slope, y0):
  """A simple linear function, applied element-wise.

  Args:
    xs (np.ndarray or float): Input(s) to the function.
    slope (float): Slope of the line.
    y0 (float): y-intercept of the line.

  """
  ys = slope * xs + y0
  return ys

セクション5: モデルについての考察 ​

チュートリアル開始からこの地点までの推定所要時間：40分

以下の質問についてグループと10分ほど議論してください：

• 以前にモデルについて聞いたことがありますか？
• これまでにモデルを作成したことはありますか？
• なぜ、Whatモデルは有用ですか？
• いつ、モデルを作成できますか？あなたの専門分野にはWhatモデルがありますか？
• モデルを構築することから学ぶことは何ですか？

フィードバックを提出する ​

まとめ ​

チュートリアルの推定所要時間：50分

このチュートリアルでは、神経データをロードし、データセットがどのように構成されているかを理解するために調査しました。次に、基本的なプロットを作成して、(1)人口全体の平均活動レベルと(2)個々のニューロンのISIの分布を視覚化しました。最後の部分では、数学的な形式を使用して生理学的な現象を理解または説明することについて考え始めました。これらの作業によって、データが「何」に似ているかを理解することができました。

これは、脳について何かを知ることができるモデルを開発するための最初のステップです。次の2つのチュートリアルで焦点を当てます。