diff --git a/.github/workflows/ebrains-push.yml b/.github/workflows/ebrains-push.yml
new file mode 100644
index 0000000000000000000000000000000000000000..a7e08d550336e76fe56c6b516755c8667f67a6b9
--- /dev/null
+++ b/.github/workflows/ebrains-push.yml
@@ -0,0 +1,31 @@
+name: Mirror to EBRAINS
+
+on:
+  push:
+    branches: [ master ]
+
+jobs:
+  sync_to_ebrains:
+    runs-on: ubuntu-latest
+    if: ${{ github.repository_owner == 'INM-6' }}
+    steps:
+      - name: Harden Runner
+        uses: step-security/harden-runner@63c24ba6bd7ba022e95695ff85de572c04a18142 # v2.7.0
+        with:
+          egress-policy: audit
+          disable-telemetry: true
+
+      - name: sycnmaster
+        uses: wei/git-sync@55c6b63b4f21607da0e9877ca9b4d11a29fc6d83
+        with:
+          source_repo: "INM-6/multi-area-model"
+          source_branch: "master"
+          destination_repo: "https://ghpusher:${{ secrets.EBRAINS_GITLAB_ACCESS_TOKEN }}@gitlab.ebrains.eu/IAS-6/multi-area-model.git"
+          destination_branch: "master"
+      - name: synctags
+        uses: wei/git-sync@55c6b63b4f21607da0e9877ca9b4d11a29fc6d83
+        with:
+          source_repo: "INM-6/multi-area-model"
+          source_branch: "refs/tags/*"
+          destination_repo: "https://ghpusher:${{ secrets.EBRAINS_GITLAB_ACCESS_TOKEN }}@gitlab.ebrains.eu/IAS-6/multi-area-model.git"
+          destination_branch: "refs/tags/*"
diff --git a/M2E_pop_rates_ref.npy b/M2E_pop_rates_ref.npy
new file mode 100644
index 0000000000000000000000000000000000000000..0dd7a7d493ff7dee24ab55f916107fb5b8566a51
Binary files /dev/null and b/M2E_pop_rates_ref.npy differ
diff --git a/M2E_statistical_test.ipynb b/M2E_statistical_test.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..705867b91f2d5f67c0ea98eebf4b867e31dc64fc
--- /dev/null
+++ b/M2E_statistical_test.ipynb
@@ -0,0 +1,168 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0bc96f39-b94f-40d5-be5c-a0d9ebdf8180",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import os\n",
+    "import sys\n",
+    "import json\n",
+    "\n",
+    "sys.path.append('./figures/MAM2EBRAINS')\n",
+    "\n",
+    "from multiarea_model import MultiAreaModel\n",
+    "from multiarea_model import Analysis\n",
+    "from M2E_compute_pop_rates import compute_pop_rates\n",
+    "\n",
+    "base_path = os.path.abspath(\".\")\n",
+    "data_path = os.path.abspath(\"simulations\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c907631c-4611-469a-941c-ddb31ef753c8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Set the threshold of matching rate of mean firing rate of all populations\n",
+    "# if smaller than this value, raise error so that the test won't pass and something is wrong and needs to be checked manually\n",
+    "match_rate_threshold = 0.95"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "97fbab67-7282-4800-a6d4-78bb3a36640a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"\n",
+    "Down-scaled model.\n",
+    "Neurons and indegrees are both scaled down to 10 %.\n",
+    "Can usually be simulated on a local machine.\n",
+    "\n",
+    "Warning: This will not yield reasonable dynamical results from the\n",
+    "network and is only meant to demonstrate the simulation workflow.\n",
+    "\"\"\"\n",
+    "d = {}\n",
+    "conn_params = {'replace_non_simulated_areas': 'het_poisson_stat',\n",
+    "               'cc_weights_factor': 1.9, # run model in Ground State\n",
+    "               'cc_weights_I_factor': 2.0}\n",
+    "network_params = {'N_scaling': 0.006,\n",
+    "                  'K_scaling': 0.006,\n",
+    "                  'fullscale_rates': os.path.join(base_path, 'tests/fullscale_rates.json')}\n",
+    "\n",
+    "sim_params = {'t_sim': 2000.,\n",
+    "              'num_processes': 1,\n",
+    "              'local_num_threads': 1}\n",
+    "\n",
+    "M = MultiAreaModel(network_params, simulation=True,\n",
+    "                   sim_spec=sim_params,\n",
+    "                   theory=True)\n",
+    "\n",
+    "M.simulation.simulate()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6730fe7b-e8bf-4aff-9497-bc614f44febd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create an instance of Analysis to load data\n",
+    "A = Analysis(M, M.simulation, data_list=['spikes'], load_areas=None)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1c907303-d2ad-4f17-83df-8b0e21749a1d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create pop_rates, load stationary firing rates, and calculate mean firing rate for all populations\n",
+    "label = M.simulation.label\n",
+    "complete_area_list = ['V1', 'V2', 'VP', 'V3', 'V3A', 'MT', 'V4t', 'V4', 'VOT', 'MSTd',\n",
+    "                      'PIP', 'PO', 'DP', 'MIP', 'MDP', 'VIP', 'LIP', 'PITv', 'PITd',\n",
+    "                      'MSTl', 'CITv', 'CITd', 'FEF', 'TF', 'AITv', 'FST', '7a', 'STPp',\n",
+    "                      'STPa', '46', 'AITd', 'TH']\n",
+    "\n",
+    "area_list = complete_area_list\n",
+    "\n",
+    "compute_pop_rates(M, data_path, label) # Compute pop_rates\n",
+    "\n",
+    "fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json') # Load stationary firing rates\n",
+    "with open(fn, 'r') as f:\n",
+    "    pop_rates = json.load(f)\n",
+    "\n",
+    "# Process spike data and calculate mean firing rate for all populations, save them in a dict\n",
+    "rates = np.zeros((len(area_list), 8))\n",
+    "for i, area in enumerate(area_list):\n",
+    "    for j, pop in enumerate(M.structure[area][::-1]):\n",
+    "        rate = pop_rates[area][pop]\n",
+    "        if rate == 0.0:\n",
+    "            rate = 1e-5\n",
+    "        if area == 'TH' and j > 3:  # To account for missing layer 4 in TH\n",
+    "            rates[i][j + 2] = rate\n",
+    "        else:\n",
+    "            rates[i][j] = rate\n",
+    "\n",
+    "rates = np.transpose(rates)\n",
+    "rates = np.array(rates)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "34d46019-38f2-4b5d-af09-575906b1d4d2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Compare the calculated mean firing rate with the reference values saved in M2E_reference_MFRs.json and calculate the percentage of matching\n",
+    "# Load the array\n",
+    "pop_rates_ref = np.load('M2E_pop_rates_ref.npy')\n",
+    "\n",
+    "# Count the percentage of cells that values match (the same)\n",
+    "comparison_result = pop_rates_ref == rates\n",
+    "# print(comparison_result)\n",
+    "same_cells = np.sum(comparison_result)\n",
+    "# print(same_cells)\n",
+    "total_cells = np.prod(pop_rates_ref.shape)\n",
+    "# print(total_cells)\n",
+    "match_rate = same_cells / total_cells\n",
+    "# print(match_rate)\n",
+    "\n",
+    "# Judge if the matching rate is lower than threshold set and raise error if yes\n",
+    "if match_rate < match_rate_threshold:\n",
+    "    raise TestError(\"Mismatch of mean firing rates over populations between the latest simulation and saved reference values, there may be problems of recent updates of dependencies, please take a check!\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "EBRAINS-23.09",
+   "language": "python",
+   "name": "ebrains-23.09"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}