Merge branch 'neurolib-dev:master' into Lead_Field_Matrix

.github/workflows/documentation.yml

@@ -29,7 +29,7 @@ jobs:
       - name: Install dependencies 🛠
         run: |
           python -m pip install --upgrade pip
-          pip install mkdocs mkdocs-material mkdocstrings mkdocstrings-python mknotebooks Pygments mkdocs-include-markdown-plugin livereload
+          pip install mkdocs mkdocs-material mkdocstrings mkdocstrings-python mknotebooks Pygments livereload
           if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
           pip install .
       - name: Build documentation 👷‍♀️

README.md

@@ -35,6 +35,8 @@ You can chose from different neural mass [models](https://github.com/neurolib-de
 
 📚 Please read the [gentle introduction](https://caglorithm.github.io/notebooks/neurolib-intro/) to `neurolib` for an overview of the basic functionality and the science behind whole-brain simulations or read the [documentation](https://neurolib-dev.github.io/) for getting started.
 
+To browse the source code of `neurolib` visit out [GitHub repository](https://github.com/neurolib-dev/neurolib).
+
 📝 <a href="#how-to-cite">Cite</a> the following paper if you use `neurolib` for your own research:
 
 > Cakan, C., Jajcay, N. & Obermayer, K. neurolib: A Simulation Framework for Whole-Brain Neural Mass Modeling. [Cogn. Comput. (2021)](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007822). 
@@ -95,6 +97,10 @@ neurolib/	 				# Main module
 ├── optimize/ 					# Optimization submodule
 	├── evolution/ 				# Evolutionary optimization
 	└── exploration/ 			# Parameter exploration
+├── control/optimal_control/			# Optimal control submodule
+	├── oc.py 				# Optimal control base class
+	├── cost_functions.py 			# cost functions for OC
+	├── /.../ 				# Implemented OC models
 ├── data/ 					# Empirical datasets (structural, functional)
 ├── utils/					# Utility belt
 	├── atlases.py				# Atlases (Region names, coordinates)
@@ -120,6 +126,7 @@ Example [IPython Notebooks](examples/) on how to use the library can be found in
 - [Example 0.6](https://mybinder.org/v2/gh/neurolib-dev/neurolib/master?filepath=examples%2Fexample-0.6-custom-model.ipynb) - Minimal example of how to implement your own model in `neurolib`
 - [Example 1.2](https://mybinder.org/v2/gh/neurolib-dev/neurolib/master?filepath=examples%2Fexample-1.2-brain-network-exploration.ipynb) - Parameter exploration of a brain network and fitting to BOLD data
 - [Example 2.0](https://mybinder.org/v2/gh/neurolib-dev/neurolib/master?filepath=examples%2Fexample-2-evolutionary-optimization-minimal.ipynb) - A simple example of the evolutionary optimization framework 
+- [Example 5.2](https://mybinder.org/v2/gh/neurolib-dev/neurolib/master?filepath=examples%2Fexample-5.2-oc-wc-model-deterministic.ipynb) - Example of optimal control of the noise-free Wilson-Cowan model
 
 A basic overview of the functionality of `neurolib` is also given in the following. 
 
@@ -282,6 +289,31 @@ This will gives us a summary of the last generation and plots a distribution of
   <img src="https://github.com/neurolib-dev/neurolib/raw/master/resources/evolution_animated.gif">
 </p>
 
+### Optimal control
+The optimal control module enables to compute efficient stimulation for your neural model. If you know how your output should look like, this module computes the optimal input. Detailes example notebooks can be found in the [example folder](https://github.com/neurolib-dev/neurolib/tree/master/examples) (examples 5.1, 5.2, 5.3, 5.4). In optimal control computations, you trade precision with respect to a target against control strength. You can determine how much each contribution affects the results, by setting weights accordingly.
+
+To compute an optimal control signal, you need to create a model (e.g., an FHN model) and define a target state (e.g., a sine curve with period 2).
+```python
+from neurolib.models.fhn import FHNModel
+model = FHNModel()
+
+duration = 10.
+model.params["duration"] = duration
+dt = model.params["dt"]
+
+period = 2.
+target = np.sin(2. * np.pi * np.arange(0, duration+dt, dt) / period)
+```
+You can then create a controlled model and run the iterative optimization to find the most efficient control input. The optimal control and the controlled model activity can be taken from the controlled model.
+```python
+model_controlled = oc_fhn.OcFhn(model, target)
+model_controlled.optimize(500) # run 500 iterations
+optimal_control = model_controlled.control
+optimal_state = model_controlled.get_xs()
+```
+
+For a comprehensive study on optimal control of the Wilson-Cowan model based on the neurolib optimal control module, see Salfenmoser, L. & Obermayer, K. Optimal control of a Wilson–Cowan model of neural population dynamics. Chaos 33, 043135 (2023). https://doi.org/10.1063/5.0144682.
+
 ## More information
 
 ### Built With
@@ -313,9 +345,11 @@ Cakan, C., Jajcay, N. & Obermayer, K. neurolib: A Simulation Framework for Whole
 
 Caglar Cakan ([email protected])  
 Department of Software Engineering and Theoretical Computer Science, Technische Universität Berlin, Germany  
-Bernstein Center for Computational Neuroscience Berlin, Germany  
+Bernstein Center for Computational Neuroscience Berlin, Germany 
 
 ### Acknowledgments
 This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) with the project number 327654276 (SFB 1315) and the Research Training Group GRK1589/2.
 
+The optimal control module was developed by Lena Salfenmoser and Martin Krück supported by the DFG project 163436311 (SFB 910).
+
 <!--end-include-in-documentation-->

docs/README.md

@@ -0,0 +1 @@
+../README.md

examples/example-0.3-fhn-minimal.ipynb

@@ -6,7 +6,7 @@
    "source": [
     "# The Fitz-Hugh Nagumo oscillator\n",
     "\n",
-    "In this notebook, the basic use of the implementation of the Fitz-Hugh Nagumo (`fhn`) model is presented. Usually, the `fhn` model is used to represent a single neuron (for example in `Cakan et al. (2014)`, \"Heterogeneous delays in neural networks\"). This is due to the difference in timescales of the two equations that define the FHN model: The first equation is often referred to as the \"fast variable\" whereas the second one is the \"slow variable\". This makes it possible to create a model with a very fast spiking mechanism but with a slow refractory period. \n",
+    "In this notebook, the basic use of the implementation of the FitzHugh-Nagumo (`fhn`) model is presented. Usually, the `fhn` model is used to represent a single neuron (for example in `Cakan et al. (2014)`, \"Heterogeneous delays in neural networks\"). This is due to the difference in timescales of the two equations that define the FHN model: The first equation is often referred to as the \"fast variable\" whereas the second one is the \"slow variable\". This makes it possible to create a model with a very fast spiking mechanism but with a slow refractory period. \n",
     "\n",
     "In our case, we are using a parameterization of the `fhn` model that is not quite as usual. Inspired by the paper by `Kostova et al. (2004)` \"FitzHugh–Nagumo revisited: Types of bifurcations, periodical forcing and stability regions by a Lyapunov functional.\", the implementation in `neurolib` produces a slowly oscillating dynamics and has the advantage to incorporate an external input term that causes a Hopf bifurcation. This means, that the model roughly approximates the behaviour of the `aln` model: For low input values, there is a low-activity fixed point, for intermediate inputs, there is an oscillatory region, and for high input values, the system is in a high-activity fixed point. Thus, it offers a simple way of exploring the dynamics of a neural mass model with these properties, such as the `aln` model.\n",
     "\n",

examples/example-0.5-kuramoto.ipynb

@@ -63,7 +63,7 @@
    "source": [
     "# Single node simulation \n",
     "\n",
-    "Here we will simulate a signal node with no noise. We then cap the phase values to be between 0 and 2*pi. We also willo plot the phase values over time."
+    "Here we will simulate a signal node with no noise. We then cap the phase values to be between 0 and 2*pi. We also will plot the phase values over time."
    ]
   },
   {
@@ -117,7 +117,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Here we simulate networks of oscillators. We will simulate a network of 8 oscillators with a global coupling strength 0.3. Here we initialize a connectivity matrix with all-to-all connectivity. We then simulate the network for 30 miliseconds assuming dt is in ms. We will also plot the phase values over time."
+    "Here we simulate networks of oscillators. We will simulate a network of 8 oscillators with a global coupling strength 0.3. Here we initialize a connectivity matrix with all-to-all connectivity. We then simulate the network for 30 milliseconds assuming dt is in ms. We will also plot the phase values over time."
    ]
   },
   {
@@ -222,13 +222,13 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Now the syncrhonization happens after 7 ms which is faster compared to the previous simulation."
+    "Now the synchronization happens after 7 ms which is faster compared to the previous simulation."
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "kuramoto",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -242,10 +242,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.16"
-  },
-  "orig_nbformat": 4
+   "version": "3.8.13"
+  }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }