The cover of my personal copy of Virchow’s Cellularpathologie, which I ordered from a second-hand book store. The book is in German, but there may be English translations available as well.

I could not resist buying a physical copy of Virchow’s Cellularpathologie. A used reproduction from 1956, nothing special. But it felt oddly appropriate. As someone still relatively new to neuroscience, I did not know this work before. However, as a side note, I had already encountered Virchow’s name almost on a daily basis: In the city where I grew up, I lived in a neighbourhood where almost all streets were named after German scientists. One of them was Virchowstraße¹ (Virchow Street), which I passed by every day on my way to school. Of course, I knew what Virchow was famous for, but I had never read any of his work.

Zoom onto the cover of Virchow’s Cellularpathologie (my personal copy).

I did not read the entire book, since it is quite long and dense, and, from a physicist’s perspective, several chapters are out of my depth. However, since it was written in German, it was of course no problem to understand the text. I found it fascinating how structured and precise this early scientific work is. It felt almost like reading a modern scientific book. Of course, the figures were not made with modern software, as we would do today. Instead, woodcut illustrations were used, which achieve such a high level of detail and clarity that they even surpass some modern illustrations in my view. The book is organized into 20 lectures, each focused on a specific aspect, and the arguments are built step by step, with careful attention to evidence and logical structure.

In this post, I want to share some of the insights I gained from engaging with Virchow’s Cellularpathologie. I will not attempt to summarize the entire book, but rather focus on some key themes and implications that I found interesting.

Virchow and the conceptual shift in medicine

Rudolf Virchow (1821–1902) is often regarded as one of the most influential figures in the history of medicine. He was a physician, pathologist, anthropologist, politician, and social reformer, and his work fundamentally changed how disease was understood in the 19th century. Today, Virchow is primarily remembered as the founder of cellular pathology, but his scientific influence extended much further. He contributed to histology, pathology, epidemiology, and public health, argued strongly for the social dimension of medicine, and introduced or popularized concepts and terms that remain central to medicine today, including thrombosis, embolism, leukemia, and neuroglia.

First pages of Virchow’s Cellularpathologie, with a photograph of Virchow himself on the left. Source: Reproduced from Rudolf Virchow, Die Cellularpathologie in ihrer Begründung auf physiologische und pathologische Gewebelehre, 1858. Photograph of the printed figure taken by myself. Original work in the public domain.

His most influential work, however, was undoubtedly Cellularpathologie, first published in 1858. The book laid the conceptual foundation for modern pathology and had profound implications for how disease itself was understood.

Virchow’s central idea is often summarized in a single statement:

Omnis cellula e cellula

Every cell arises from another cell.

At first glance, this sounds almost trivial from a modern perspective. However, in the mid 19th century, this statement directly contradicted older ideas about spontaneous generation and the formation of tissues from amorphous substances. It established continuity at the cellular level as the fundamental principle of life. More importantly, Virchow extended this principle from physiology to pathology. Disease, in his view, is not a disturbance of abstract bodily “fluids” or humors such as blood, phlegm, black bile, and yellow bile, as in ancient medicine, but a disturbance of cells. This is the conceptual core of cellular pathology.

This insight has profound implications: If every physiological process is rooted in cellular processes, then every pathological process must also be reducible to changes in cellular structure, function, or proliferation. In modern terms, this marks a shift from descriptive medicine at the level of organs and symptoms toward a mechanistic and microscopic understanding of disease.

The lectures

Virchow’s Cellularpathologie is not a monolithic treatise in the traditional sense. It is a collection of 20 lectures, delivered in early 1858 at the Pathological Institute in Berlin and later transcribed and edited for publication. This also explains the style of the work: It is didactic, structured around demonstrations, and closely tied to visual material such as microscopic preparations and woodcut illustrations.

Cover of the first lecture of Virchow’s Cellularpathologie. Source: Woodcut illustration reproduced from Rudolf Virchow, Die Cellularpathologie in ihrer Begründung auf physiologische und pathologische Gewebelehre, 1858. Photograph of the printed figure taken by myself. Original work in the public domain.

The work begins with a systematic overview of cell theory, integrating botanical and zoological observations. The inclusion of plant cells alongside animal tissues reflects the universality of the cellular principle. The same structural logic applies across living systems, whether in plant tissue, cartilage, liver, or nervous tissue, as illustrated in the early figures of the book.

Woodcut illustration of an injection preparation of the skin, vertical cut. E: Epidermis, R: Rete Malpighii, P: Papilla with its up and down projecting vessels, C: Cutis. Source: Woodcut illustration (Fig. 44) reproduced from Rudolf Virchow, Die Cellularpathologie in ihrer Begründung auf physiologische und pathologische Gewebelehre, 1858. Photograph of the printed figure taken by myself. Original work in the public domain.

From this foundation, Virchow gradually develops the framework of cellular pathology in detail. Much of the book can be understood as a systematic attempt to establish the cell as the central explanatory unit of medicine.

Woodcut illustration of Ganglion cells from the central organs: A, B, C from the spinal cord, based on preparations by Mr. Gerlach; D from the cerebral cortex. A: Large, multi-dendritic (multipolar, polyclonal) cells from the anterior horns (motor cells). B: Smaller cells with 3 larger processes from the posterior horns (sensory cells). C: Two-branched (bipolar, diclonal), more rounded cells from the vicinity of the posterior commissure (sympathetic cells). Magnification: 300. Source: Woodcut illustration (Fig. 89) reproduced from Rudolf Virchow, Die Cellularpathologie in ihrer Begründung auf physiologische und pathologische Gewebelehre, 1858. Photograph of the printed figure taken by myself. Original work in the public domain.

In the lectures, Virchow explicitly argues against older humoral and organ-centered views of disease and instead localizes pathological processes at the cellular level. However, what is particularly interesting here is how concretely he develops this argument. Rather than presenting a purely abstract theory, he continuously moves between microscopic observations, tissue morphology, pathological examples, and conceptual interpretation.

Cover of the 12th lecture of Virchow’s Cellularpathologie, which is focused on the nervous system. Source: Woodcut illustration reproduced from Rudolf Virchow, Die Cellularpathologie in ihrer Begründung auf physiologische und pathologische Gewebelehre, 1858. Photograph of the printed figure taken by myself. Original work in the public domain.

Similarly, Virchow repeatedly emphasizes that pathological processes are not fundamentally separate from physiological ones. Disease is not treated as something metaphysically distinct from normal life processes, but as an alteration, dysregulation, or exaggeration of normal cellular activity. In his lectures, he establishes this view not only philosophically, but through numerous examples involving growth, degeneration, inflammation, and tissue remodeling. The detailed discussion of inflammation is especially interesting here: Inflammation is no longer described merely through externally visible symptoms such as heat, redness, and swelling, but as a sequence of cellular events involving the accumulation, transformation, and activity of specific cells. Likewise, tumors are interpreted not as foreign entities or parasitic growths, but as the result of abnormal cellular proliferation originating from the body’s own tissues.

And Virchow reinforces his arguments with an, at that time, extensive use of microscopic images, which range from plant cells to cartilage, liver cells, capillaries, ganglion cells, and pathological tissues:

Schematic representation of nerve fiber patterns in the cerebellar cortex according to Gerlach. (Microscopic Studies, Plate I, Fig. 3.) A: white matter; B, C: gray matter; B: granular layer; C: cell layer. Source: Woodcut illustration (Fig. 91) reproduced from Rudolf Virchow, Die Cellularpathologie in ihrer Begründung auf physiologische und pathologische Gewebelehre, 1858. Photograph of the printed figure taken by myself. Original work in the public domain.

Another important aspect is Virchow’s insistence on continuity. There is no abrupt transition from health to disease at the level of fundamental principles. Instead, there exists a continuous spectrum of cellular states and transformations. This idea feels remarkably modern, since many diseases today are understood not as entirely separate biological phenomena, but as dysregulated versions of normal cellular processes.

Finally, and in direct connection to our earlier discussion of Helmholtz, Virchow introduced the concept of neuroglia as a distinct cellular component of nervous tissue. This distinction between neurons and glial cells provided an important conceptual framework for later neuroscience and fundamentally changed how the nervous system was understood at the cellular level. What Virchow recognized was that the nervous system could not be reduced to nerve fibers and nerve cells alone. It also contained another cellular and structural component, which he interpreted as a kind of connective or supporting tissue of the nervous system. The term neuroglia itself expresses this idea: A “nerve glue”, a tissue that holds the nervous elements together:

Woodcut illustration of the ependyma of the ventricles and neuroglia from the floor of the fourth ventricle. E: epithelium, N: nerve fibers. In between lies the free portion of the neuroglia, containing numerous connective tissue cells and nuclei; at v is a blood vessel; elsewhere, numerous corpora amylacea are visible, which are shown in isolation at ca. Magnification: 300. Source: Woodcut illustration (Fig. 94) reproduced from Rudolf Virchow, Die Cellularpathologie in ihrer Begründung auf physiologische und pathologische Gewebelehre, 1858. Photograph of the printed figure taken by myself. Original work in the public domain.

From a modern perspective, this interpretation was still limited, because glial cells are no longer understood as merely passive support elements. Astrocytes, oligodendrocytes, microglia, and other glial cell types are now known to participate actively in homeostasis, metabolism, myelination, immune defense, synaptic regulation, and disease processes. But Virchow’s conceptual step was nevertheless crucial: He separated the nervous system into more than one cellular category and thereby opened the possibility of asking what these non-neuronal cells are, where they come from, and what they do.

Conceptual implications

What I take away from researching and reading Virchow’s Cellularpathologie is that it represents a major conceptual shift in the history of medicine. It is not just a collection of observations or a description of pathological phenomena. It is an attempt to reorganize the entire field of medicine around the cell as the fundamental unit of life and disease. Virchow constructs a coherent explanatory framework that connects anatomy, physiology, pathology, microscopy, and tissue structure into a unified view of disease. In this sense, the book can be regarded as on of the scientifically most important works in the history of medicine, and it laid the groundwork for many subsequent developments in pathology, histology, and neuroscience.

Several aspects of Virchow’s reasoning still feel surprisingly familiar, at least from what I can recognize as a non-medical reader. The idea that disease emerges from altered cellular states, that pathological and physiological processes are deeply connected, and that microscopic structures can provide mechanistic explanations for macroscopic symptoms remains central to biomedical science today.

The garden pavilion of the Juliusspital in Würzburg housed Virchow’s Institute of Pathology until 1853. Source: Wikimedia Commonsꜛ (license: CC BY-SA 4.0).

Engaging with Virchow after reading Helmholtz’s dissertation also highlights an interesting historical contrast. Helmholtz’s dissertation is to some extent a bit more observational and comparative. Virchow, by contrast, is much more focused on conceptual integration and theoretical organization. He builds a framework in which such structures acquire general pathological meaning. In this sense, the two works together provide an interesting glimpse into the emergence of modern life sciences during the 19th century, which could, in a very simplified way, be seen as: One focused on observation and comparative anatomy, the other on conceptual unification and medical interpretation.

This little excursion into the history of science was quite refreshing. Ideas we now take for granted, not only in neuroscience and medicine, must have been revolutionary at that time and had to be established through careful argumentation and evidence, which is not so different from today. I also enjoyed reading about the biographies of 19th and early 20th century² scientists, who were often polymaths³ and engaged in multiple fields of science, medicine, and even politics. This again reminds me that the boundaries between disciplines were much more fluid at that time, and that scientific progress often involved contributions from individuals with broad interests and expertise. An attitude that I think is worth remembering today.

References and further reading

• Rudolf Virchow, Die Cellularpathologie in ihrer Begründung auf physiologische und pathologische Gewebelehre: Zwanzig Vorlesungen, gehalten während der Monate Februar, März und April 1858 im pathologischen Institute zu Berlin, 1858, Verlag von August Hirschwald, Berlin; online available for free at MDZ Digitale Sammlungenꜛ or on archive.orgꜛ (original work in the public domain)
• Julia Heideklang, H.-J. Pflüger, Helmut Kettenmann, De fabrica systematis nervosi evertebratorum. Die kommentierte Dissertation von Hermann Helmholtz, 2021, Wbg Academic, ISBN: 9783534400942, online PDFꜛ, Websiteꜛ

Footnotes

My personal hard copy of De fabrica systematis nervosi evertebratorum. Die kommentierte Dissertation von Hermann Helmholtz. The book contains the original Latin text, a German translation, an English translation, and extensive commentary by Julia Heideklang, Hans-Joachim Pflüger, and Helmut Kettenmann. The fact that Helmholtz’s dissertation dealt with the nervous system was quite surprising to me, since I had always associated him with physics and mathematics.

And if this was not surprising enough, Kettenmann made a further remark that was even more striking: Helmholtz’s dissertation was an anatomical study of invertebrate nervous systems, written at a time when the basic conceptual distinction between neurons and glial cells did not yet exist. In other words, Helmholtz investigated the structure of the nervous system without the conceptual framework that we now take for granted. This is remarkable, since the distinction between neurons and glia was introduced only later by Rudolf Virchow in 1856/1858. Helmholtz was therefore working in a conceptual landscape that was still largely undefined, and yet he was already engaging with the cellular organization of nervous tissue.

Thanks to the work of Helmut Kettenmann and his two colleagues, this early scientific text is now accessible again. The translation is freely available onlineꜛ, and I could not resist buying a hard copy of the book to hold in my hands and read through it. In this post, I would therefore like to share some impressions and insights from reading the dissertation, and to reflect on its significance, not as a historical curiosity, but as an early contribution to neuroscience that has largely been forgotten.

Helmholtz: Biography and what we usually know about him

Hermann von Helmholtz (1821–1894) is primarily known as one of the central figures of 19th century physics. For physicists, his name is associated with a range of foundational results across electrodynamics, thermodynamics, and mathematical physics.

Left: Hermann von Helmholtz in 1848. Source: Wikimedia Commonsꜛ (license: public domain). Right: Last photograph of von Helmholtz, taken three days before his final illness. Source: Wikimedia Commonsꜛ (license: CC0 1.0 Universal Public Domain Dedication). Helmholtz is widely regarded as one of the greatest physicists of the 19th century, and his contributions to physics are foundational. However, his early work on the nervous system is not commonly known and is often overlooked in standard accounts of his scientific career. This is partly because his dissertation was written in Latin and therefore remained largely inaccessible for a long time, but it also reflects a broader tendency in the history of science to reduce scientists to only one part of their intellectual work in my view.

One of the most familiar formulations associated with Helmholtz is the Helmholtz decomposition theorem (in German often just referred to as the Helmholtz theorem). Let’s briefly recall what this is about. In vector calculus, the theorem states that any sufficiently smooth vector field can be decomposed into a divergence-free (solenoidal) part and a curl-free (irrotational) part. Mathematically, for a vector field $\mathbf{F}(\mathbf{x})$ that vanishes at infinity, we can write:

Here, $\phi$ is a scalar potential and $\mathbf{A}$ a vector potential. Physically, this means that any field can be separated into a component that originates from sources or sinks and a component that describes rotational or circulating structures. This decomposition underlies much of classical electrodynamics and fluid dynamics, where electric, magnetic, or velocity fields are naturally described in terms of potentials.

Closely related is the Helmholtz equation, which appears throughout wave physics:

Here, $\psi$ describes a spatial wave field, such as a sound pressure distribution, an electromagnetic field amplitude, or a quantum mechanical wavefunction. The Laplacian $\nabla^2$ captures how the field changes spatially, while the parameter $k$ is the wave number, related to the wavelength $\lambda$ by

The Helmholtz equation emerges when oscillatory systems are studied in the frequency domain. Instead of describing how a wave evolves in time, it describes the spatial structure of standing or propagating waves at a fixed frequency. Solutions therefore represent spatial wave patterns and resonant modes. The equation appears in acoustics, optics, electrodynamics, and quantum mechanics, but also in plasma physics and space physics, where wave-like disturbances propagate through magnetized plasmas. Plasma waves in the solar wind, magnetospheric oscillations, and electromagnetic wave propagation in ionized media can all lead, under suitable approximations, to Helmholtz-type equations.

Left: Two sources of radiation in the plane, represented mathematically by a source function $f$, which vanishes in the blue region. Right: The real part of the resulting field $A$. Here, $A$ is the solution of the inhomogeneous Helmholtz equation $(\nabla^2 + k^2)A = -f$. Source: Wikimedia Commonsꜛ and Wikimedia Commonsꜛ (license: public domain). The Helmholtz equation describes how localized sources generate spatial wave fields and interference patterns. Such equations form part of the mathematical foundation of modern wave physics.

Helmholtz’s scientific work, however, extended far beyond vector calculus and wave theory. He also formulated the principle of conservation of energy in a modern quantitative framework, thereby contributing to the unification of mechanics, heat theory, and electromagnetism. At the same time, he remained deeply connected to physiology and medicine throughout his career. In experimental physiology, for example, he famously measured the finite speed of nerve conduction, thereby refuting the older assumption that neural signaling occurred instantaneously. This finding is already remarkable for itself, as at his time, the experiments by Hodgkin and Huxley that would later provide a detailed mechanistic explanation of nerve conduction were still decades away. But it also illustrates how Helmholtz’s scientific interests spanned both physics and physiology, and how he applied rigorous quantitative methods to the study of biological systems.

So, Helmholtz did not simply move from medicine into physics and leave biology behind. Rather, he approached biological systems with the same quantitative and mechanistic mindset that characterized his physical work. In many ways, he helped establish the broader 19th century idea that physiological processes could be investigated with the tools and methods of physics.

From a physicist’s perspective, Helmholtz is therefore a unifying figure who connects field theory, wave physics, thermodynamics, and physiology. The idea that this same person began his scientific career with a detailed anatomical study of invertebrate nervous systems is therefore unexpected at first glance, but perhaps also strangely fitting.

Helmholtz’s dissertation

Helmholtz’s 1842 dissertation, De fabrica systematis nervosi evertebratorum, is situated in a specific scientific context. At that time, microscopy had only recently reached a level that allowed the visualization of cellular structures. Christian Ehrenberg (1795-1876)ꜛ had published the first images of neuronsꜛ (Ehrenberg (1832, 1836); Chvátal (2015)ꜛ), including ganglia of the leech, and the emerging cell theory of Schleidenꜛ and Schwannꜛ proposed that all tissues are composed of cells in 1838/1839 (Schwann (1838); Schwann (1839)).

First image of a neuron: A drawing of the leech ganglion by Christian Ehrenberg (1836). Individual magnified neurons are displayed on the left. Source: Julia Heideklang, H.-J. Pflüger, Helmut Kettenmann, De fabrica systematis nervosi evertebratorum. Die kommentierte Dissertation von Hermann Helmholtz, 2021, page 10, (license: CC BY-SA 4.0).

Top Left: Theodore Schwann (1810-1882); from The Story of Nineteenth-Century Science, page 377, 1904, by Henry Smith Williams. Source: Wikimedia Commonsꜛ (license: public domain). Top Right: Matthias Schleiden (1804-1881), between 1882 and 1883. Source: Wikimedia Commonsꜛ (license: public domain). Bottom: Human cancer cells with nuclei (specifically the DNA) stained blue. The central and rightmost cell are in interphase, so the entire nuclei are labeled. The cell on the left is going through mitosis and its DNA has condensed. Source: Wikimedia Commonsꜛ (license: public domain). – Schwann and Schleiden are credited with the formulation of the cell theory, which posits that all living organisms are composed of cells. This was a foundational concept in biology, and it set the stage for later investigations into the cellular structure of tissues, including the nervous system. Helmholtz’s dissertation can be seen as an early attempt to apply these emerging ideas to the study of nervous tissue in invertebrates.

However, while vertebrate nervous systems had begun to be described, a systematic investigation of invertebrate nervous systems was largely missing. Helmholtz’s supervisor, Johannes Müllerꜛ, assigned precisely this problem to the young medical student: To determine whether the structural principles of nervous systems are conserved across species.

Experimental animals used by Helmholtz in his dissertation: Giant slug (Arion empiricorum; today: Arion ater), Great pond snail (Limnaea stagnalis), Burgundy snail (Helix pomatia), Great ramshorn snail (Planorbis corneus; today: Planorbarius corneus), Freshwater pearl mussel (Unio margaritifera; today: Margaritifera margaritifera), Common leech (Hirudo vulgaris), Earthworm (Lumbricus terrestris), Medicinal leech (Hirudo medicinalis). Source: Julia Heideklang, H.-J. Pflüger, Helmut Kettenmann, De fabrica systematis nervosi evertebratorum. Die kommentierte Dissertation von Hermann Helmholtz, 2021, page 14, (license: CC BY-SA 4.0).

Helmholtz approached this question through comparative anatomical analysis of a range of invertebrates, including leeches, worms, crustaceans, and insects. His method was fundamentally microscopic dissection and observation. For example, he describes how a nerve or nerve cord, taken from a living animal and mechanically separated on a glass plate, reveals individual nerve fibers under appropriate preparation. Similarly, by dissecting ganglia, he identifies cellular structures, including what he describes as “caudate cells”, that is, cell bodies with extensions.

Nervous system of the Medicinal leech. (A) Scheme of the nervous system; (B) Upper part of the supraesophageal ganglion (brain); (C) Fourth ganglion of the ventral nerve cord; Retzius 1891. Helmholtz describes the nervous system of the leech in detail, identifying nerve cords, ganglia, and cellular structures. This figure from Retzius (1891) illustrates the organization of the leech nervous system, which Helmholtz would have studied in his dissertation. Source: Julia Heideklang, H.-J. Pflüger, Helmut Kettenmann, De fabrica systematis nervosi evertebratorum. Die kommentierte Dissertation von Hermann Helmholtz, 2021, page 19, (license: CC BY-SA 4.0).

Thus, Helmholtz is effectively describing neurons as entities consisting of a cell body and processes, without yet having the conceptual framework of the neuron doctrine (the principle that the nervous system consists of distinct individual neurons rather than a continuous network). The terminology is not modern, but the observational content clearly anticipates it.

A central structural motif in his analysis is the organization of nerve cords and ganglia. He identifies longitudinal nerve cords connected by ganglia and observes that these ganglia contain lobed structures receiving nerve fibers from the cords. He further describes plexus-like arrangements, where multiple nerve branches interconnect and exchange fibers, in a manner analogous to plexuses in vertebrates.

One of the most significant conceptual conclusions of the dissertation is the assertion of structural homology between invertebrate and vertebrate nervous systems. Helmholtz argues that, despite differences in overall organization, the fundamental elements are the same. Nerve fibers, cellular bodies, and their interconnections follow a common structural plan across species. This is explicitly reflected in his comparison of branching patterns in invertebrate nerve plexuses to those found in vertebrates, where each branch carries fibers corresponding to those entering and leaving connected structures.

Nervous system of the crayfish after Huxley (1880): (A) Scheme of the nervous system (Retzius 1890), (B) Supraesophageal ganglion (brain), (Retzius 1890), (C) Third thoracic ganglion (Retzius 1890), (D) Individual nerve cells (Retzius 1890), (E) First abdominal ganglion (Retzius 1890). Helmholtz argues in his dissertation that the structural organization of the nervous system in invertebrates, such as the crayfish, shares fundamental features with vertebrate nervous systems: Nerve fibers, ganglia, and plexus-like arrangements are common motifs that reflect a shared structural plan. Source: Julia Heideklang, H.-J. Pflüger, Helmut Kettenmann, De fabrica systematis nervosi evertebratorum. Die kommentierte Dissertation von Hermann Helmholtz, 2021, page 21, (license: CC BY-SA 4.0).

In modern terms, this can be interpreted as an early formulation of a unifying principle of nervous system organization: That neural circuits are constructed from recurring cellular and connective motifs, independent of the organism. This is not yet an evolutionary argument in the Darwinian sense, but it is clearly a structural generalization across species.

Another important aspect is actually methodological. Helmholtz relies on careful mechanical dissection, microscopy, and comparative reasoning. There is no electrophysiology, no staining techniques in the modern sense, and no formal theory of synaptic connectivity. Yet, within these constraints, he identifies consistent structural patterns and interprets them in a general framework.

It is also crucial to emphasize what is absent. As mentioned in the introduction, there is no distinction between neurons and glial cells, because this conceptual separation had not yet been introduced. The nervous system is treated as a collection of fibers and cellular elements without a clear functional differentiation between support and signaling components. This absence highlights how early Helmholtz’s work is situated within the development of neuroscience.

What also surprised me is the absence of any figures in the dissertation. The text in this new translation is entirely descriptive, both in the original Latin version and in the translation, and while it contains detailed descriptions of structures, there are no illustrations or diagrams. This is quite different from what we know from today’s scientific publications, where figures play a central role in conveying information. At the moment, I can’t rule out whether this was an editorial choice in the translation or whether the original dissertation also lacked figures. But either way, it emphasizes the textual and descriptive nature of scientific communication at that time. If you have already seen other scientific texts from the early 19th century, you might know that this was not uncommon. The detailed descriptions of concepts and observations were conveyed entirely through words, which is a significant difference from modern scientific writing.

Taken together, the dissertation can be understood as an early comparative neuroanatomical study that establishes three key points:

1. First, nervous systems across different animal classes share a common structural basis.
2. Second, these systems are composed of cellular elements and fibers that can be identified microscopically.
3. Third, their organization into ganglia, cords, and plexuses follows systematic and reproducible patterns.

These insights, while not yet framed in the language of modern neuroscience, clearly anticipate later developments in the understanding of nervous system structure and organization. Helmholtz’s work can therefore be seen as a foundational contribution to the field, even if it has been largely overlooked in standard historical accounts.

Conclusion

What remains, after reading this dissertation in light of Helmholtz’s later work, is a certain sense of discontinuity that is, in fact, highly revealing. The figure who would later formalize energy conservation, develop mathematical frameworks for fields and waves, and shape theoretical physics, begins with meticulous anatomical observations of leeches and worms under a microscope.

From a modern perspective, this is not merely an anecdotal curiosity. It reflects a broader feature of 19th century science (and even in earlier times), where disciplinary boundaries were not yet rigidly established. The same individual could move from cell structure to electrodynamics without the conceptual friction that would exist today. Other famous figures such as Carl Friedrich Gauss, Michael Faraday, and James Clerk Maxwell also had similarly broad scientific interests that spanned multiple fields.

At the same time, the dissertation itself is more than a historical footnote. It contains a clear recognition of structural unity across nervous systems and an implicit cellular perspective that anticipates later developments. That this work was written in Latin and remained largely inaccessible for nearly two centuries may explain why it has not entered the standard narrative of neuroscience.

And Kettenmann’s remark on this work, in retrospect, points to something more general: The early history of neuroscience was less linear than our retrospective narratives often suggest. Ideas that now seem closely connected, such as cell theory, comparative neuroanatomy, glia, neurons, and physiological mechanism, did not emerge as one coherent program. They developed in fragments, in different disciplines, with different terminologies and different methods. Helmholtz’s early contribution to neuroscience was neither trivial nor naive. It was one such fragment. Revisiting it does not change the established milestones of neuroscience, but it does refine the picture. It shows that some of the core ideas about the cellular organization and comparative structure of nervous systems were already being articulated at a time when even the basic terminology had not yet been established.

Personally, this is also why I find Helmholtz’s scientific trajectory so interesting. I began in physics myself and later moved into neuroscience, although, of course, in a completely different historical and scientific context. For that reason, the trajectory from physical thinking to biological questions feels familiar, and reminds me that the border between physics and neuroscience is not as sharp as it sometimes appears. Helmholtz’s dissertation makes this visible in an unusually early and concrete form: The nervous system was already a place where anatomy, physiology, microscopy, and physical reasoning could meet.

As mentioned at the beginning, the new translation is sharedꜛ under a Creative Commons license and is freely available onlineꜛ. For anyone interested in the history of neuroscience, physiology, or 19th century science more broadly, it is well worth reading. Feel free to share your thoughts and insights in the comments below.

References and further reading

• Julia Heideklang, H.-J. Pflüger, Helmut Kettenmann, De fabrica systematis nervosi evertebratorum. Die kommentierte Dissertation von Hermann Helmholtz, 2021, Wbg Academic, ISBN: 9783534400942, online PDFꜛ, Websiteꜛ
• Chvátal, A., Discovering the Structure of Nerve Tissue: Part 1: From Marcello Malpighi to Christian Berres, 2015, Journal of the History of the Neurosciences, 24(3), 268–291. doi: 10.1080/0964704X.2014.977676ꜛ
• Ehrenberg, C. G., Über das neueste Mikroskop, von Pistor und Schiek in Berlin, gefertigt im Januar 1832, 1832, Ein Schreiben von Hrn. Ehrenberg an den Herausgeber, Januarheft 1832 der Analen der Physik, 188-192
• Ehrenberg, C. G., Abhandlung der Akademie: Beobachtung einer bisher unbekannren auffallenden Structur des Seelenorgans bei Menschen und Thieren, 1836, Berlin
• Schleiden, M., Beträge zur Phytogenesis, 1838, Archiv für Anatomie, Physiologie und wissenschaftliche Medicin, Weit et Comp. Berlin
• Schwann, T., Mikroskopische Untersuchungen über die Übereinstimmung in der Struktur und dem Wachstum der Thiere und Pflanzen, 1839, Sander’sche Buchhandlung, Berlin
• Rudolf Virchow, Die Cellularpathologie in ihrer Begründung auf physiologische und pathologische Gewebelehre: Zwanzig Vorlesungen, gehalten während der Monate Februar, März und April 1858 im pathologischen Institute zu Berlin, 1856, Verlag von August Hirschwald, Berlin

Result of a spike pairing simulation using the Clopath rule. The plot shows the change in synaptic weight as a function of the relative timing of pre- and postsynaptic spikes. The Clopath rule captures both potentiation and depression of synaptic weights based on the timing of spikes and the postsynaptic membrane potential. We will simulate this experiment in the section “Spike pairing simulation” below.

This makes the Clopath rule especially important in computational neuroscience. It provides a compact phenomenological framework that links synaptic plasticity to membrane dynamics, spike timing, and homeostatic stabilization. As a result, it is widely used in spiking neural network simulations that aim to remain closer to biological plasticity mechanisms than pair-based STDP rules.

Mathematical description

The central idea of the Clopath rule is that potentiation and depression depend on different combinations of presynaptic activity and postsynaptic voltage. In contrast to simple STDP rules, the postsynaptic membrane potential is not treated as an implicit byproduct of spike times, but as a relevant dynamical quantity in its own right.

A simplified schematic form of the rule can be written as

where

• $w_{ij}$ is the synaptic weight from presynaptic neuron $i$ to postsynaptic neuron $j$,
• $A_\text{LTD}$ and $A_\text{LTP}$ set the relative strength of the depressive and potentiating plasticity components,
• $x_i(t)$ and $\bar{x}_i(t)$ are presynaptic traces,
• $u_j(t)$ is the postsynaptic membrane potential,
• $\bar{u}_j(t)$ is a low pass filtered postsynaptic voltage,
• $\theta_-$ and $\theta_+$ are voltage thresholds for depression and potentiation, respectively, and
• $[\,\cdot\,]_+ = \max(\cdot,0)$ denotes rectification.

This expression captures the main logic of the model:

• depression is triggered when presynaptic activity coincides with sufficient postsynaptic depolarization above a lower threshold
• potentiation requires stronger postsynaptic depolarization above a higher threshold and depends on presynaptic activation through a trace term
• the two processes are therefore asymmetric in their voltage dependence

The precise formulation differs somewhat between the original papers and specific simulator implementations such as NEST (see below), but the conceptual structure remains the same.

A typical presynaptic trace is described by

where $\tau_x$ is the decay time constant and $t_i^f$ denotes the times of presynaptic spikes. $\delta(t - t_i^f)$ represents the Dirac delta function indicating a presynaptic spike at time $t_i^f$. $f$ denotes the time of the presynaptic spike. The (low-pass) filtered postsynaptic voltage potential can be written as

where $\tau_u$ is the corresponding voltage filter (decay) time constant.

In addition to its voltage dependence, the original Clopath model also contains a homeostatic component that prevents runaway synaptic growth. This is one of the key reasons why the rule is more stable in recurrent network simulations than many simpler STDP variants.

The synaptic weights are typically bounded such that

In some versions of the model, the LTD amplitude is further modulated by a homeostatic postsynaptic activity term, for example in the schematic form

where $\alpha$ is a scaling factor, $u_\text{ref}$ is a reference level, and $\bar{u}_{j,\mathrm{homeo}}$ denotes a slowly varying postsynaptic activity or voltage measure rather than the instantaneous membrane potential itself. The exact form depends on the specific model variant or simulator implementation. This allows the rule to counteract runaway potentiation and helps maintain stable learning in recurrent networks.

Implementation in neural network simulations

To use the Clopath rule in a spiking network simulation, one must supply not only pre- and postsynaptic spikes but also the postsynaptic membrane voltage and its filtered versions. For this reason, the rule is often coupled to a neuron model that provides a realistic subthreshold voltage trajectory. A common choice is the adaptive exponential integrate and fire (AdEx) neuron model, which captures leak dynamics, exponential spike initiation, and adaptation. A standard AdEx membrane equation is

where $u(t)$ is the membrane potential, $C$ is the membrane capacitance, $g_L$ is the leak conductance, $E_L$ is the resting potential, $\Delta_T$ controls the sharpness of spike initiation, $V_T$ is the (adaptive) threshold parameter, $I_\text{adapt}(t)$ is the adaptation current, $I_\text{syn}(t)$ is the synaptic current, and $I_\text{ext}(t)$ is the external input current.

The adaptation current is commonly modeled as

where $a$ is the adaptation conductance and $\tau_\text{adapt}$ is the adaptation time constant.

In the context of the Clopath rule, the synaptic current $I_\text{syn}(t)$ is given by the sum of inputs from all presynaptic neurons, weighted by their synaptic weights $w_{ij}$:

In this setting, the Clopath synapse does not replace the neuron model. Rather, it uses the neuron model’s voltage trajectory to determine when synaptic potentiation or depression should occur. This separation is important conceptually: The AdEx model describes membrane dynamics, whereas the Clopath rule describes plasticity at synapses.

In some AdEx-based implementations, the threshold potential is itself adaptive and evolves dynamically after each spike. A simple phenomenological form is

where $V_{T,\infty}$ is the baseline threshold value and $\tau_{V_T}$ is the corresponding relaxation time constant.

In Clopath and Gerstner (2010)ꜛ, the interaction between spike timing and voltage was further analyzed in a formulation that includes a depolarizing spike afterpotential. This can be represented by an additional variable $z(t)$ that is incremented at spike times and otherwise decays exponentially:

Such additions are useful because they help reproduce experimental pairing protocols more faithfully.

Why the Clopath rule matters

The Clopath learning rule is significant for several reasons:

Biological plausibility

A major advantage of the Clopath rule is that it links synaptic modification to the postsynaptic voltage state. This reflects the fact that in biological neurons, synaptic plasticity is shaped not only by spike timing but also by dendritic and somatic depolarization. In this sense, the rule is substantially closer to experimentally observed plasticity than purely timing-based STDP models.

Voltage dependence

The model distinguishes between weaker depolarization regimes associated with depression and stronger depolarization regimes associated with potentiation. This captures an important principle of real synapses: The effect of presynaptic input depends on the electrical state of the postsynaptic cell. The same presynaptic spike train can therefore induce very different weight changes depending on whether the postsynaptic neuron is weakly depolarized, strongly depolarized, or near firing threshold.

Homeostasis

The Clopath model includes homeostatic ingredients that help constrain weight growth and maintain stable activity over time. This is especially important in recurrent networks (RNN), where unconstrained Hebbian plasticity can otherwise produce runaway excitation or pathological synchrony.

Computational usefulness

From a modeling perspective, the Clopath rule is attractive because it remains relatively compact while capturing several important experimental dependencies. In particular, it provides a practical compromise between biological realism and computational tractability: It goes beyond purely timing-based STDP rules by incorporating postsynaptic voltage, yet remains simple enough to be used in larger spiking network simulations. This makes it useful in studies of receptive field formation, self-organization, synaptic competition, and activity-dependent structure in recurrent circuits. Because each synapse is updated locally on the basis of presynaptic activity and the postsynaptic voltage state, even reciprocal connections between two neurons do not have to evolve symmetrically. In this way, the rule can contribute to the emergence of structured and direction-specific connectivity patterns rather than merely uniform global weight changes. Whenever a plasticity model is needed that depends on more than spike timing alone, the Clopath rule is a natural candidate.

Spike pairing simulation

The Clopath rule (clopath_synapseꜛ) and the corresponding voltage-based neuron model (aeif_psc_delta_clopathꜛ) are both implemented in the NEST simulator. To illustrate the behavior of the rule, we can reproduce the NEST tutorial “Clopath Rule: Spike pairing experiment”ꜛ. The goal is to study how synaptic weights change under repeated spike pairings with different temporal structure and different effective pairing frequencies.

Let’s start by importing the necessary libraries:

import matplotlib.pyplot as plt
import numpy as np
import nest

# set the verbosity of the NEST simulator:
nest.set_verbosity("M_WARNING")

# Set global properties for all plots
plt.rcParams.update({'font.size': 12})
plt.rcParams["axes.spines.top"]    = False
plt.rcParams["axes.spines.bottom"] = False
plt.rcParams["axes.spines.left"]   = False
plt.rcParams["axes.spines.right"]  = False

nest.ResetKernel()

# define the simulation resolution:
resolution = 0.1

Next, we define the parameters of the neuron.,

# define the parameters of the neuron:
clopath_neuron_params = {
    "V_m": -70.6,           # [mV] membrane potential
    "E_L": -70.6,           # [mV] leak reversal potential
    "C_m": 281.0,           # [pF] membrane capacitance
    "theta_minus": -70.6,   # [mV] threshold for u (LTD)
    "theta_plus": -45.3,    # [mV] threshold for u_bar (LTP)
    "A_LTD": 14.0e-5,       # [nS] amplitude of LTD
    "A_LTP": 8.0e-5,        # [nS] amplitude of LTP
    "tau_u_bar_minus": 10.0,# [ms] time constant for u_bar_minus (LTD)
    "tau_u_bar_plus": 7.0,  # [ms] time constant for u_bar_plus (LTP)
    "delay_u_bars": 4.0,    # [ms] delay with which u_bar_[plus/minus] are processed to compute the synaptic weights.
    "a": 4.0,               # [nS] subthreshold adaptation
    "b": 0.0805,            # [pA] spike-triggered adaptation
    "V_reset": -70.6 + 21.0,# [mV] reset potential
    "V_clamp": 33.0,        # [mV] clamping potential
    "t_clamp": 2.0,         # [ms] clamping time
    "t_ref": 0.0,           # [ms] refractory period
}

We then need to define the spike times for the pre- and postsynaptic neurons. They are arranged in two lists, spike_times_pre and spike_times_post, each containing the spike times for different pairs of pre- and postsynaptic neurons. The simulation will run for each pair of spike times:

spike_times_pre = [
    # Presynaptic spike before the postsynaptic
    [20.0, 120.0, 220.0, 320.0, 420.0],
    [20.0, 70.0, 120.0, 170.0, 220.0],
    [20.0, 53.3, 86.7, 120.0, 153.3],
    [20.0, 45.0, 70.0, 95.0, 120.0],
    [20.0, 40.0, 60.0, 80.0, 100.0],
    # Presynaptic spike after the postsynaptic
    [120.0, 220.0, 320.0, 420.0, 520.0, 620.0],
    [70.0, 120.0, 170.0, 220.0, 270.0, 320.0],
    [53.3, 86.6, 120.0, 153.3, 186.6, 220.0],
    [45.0, 70.0, 95.0, 120.0, 145.0, 170.0],
    [40.0, 60.0, 80.0, 100.0, 120.0, 140.0]]

spike_times_post = [
    [10.0, 110.0, 210.0, 310.0, 410.0],
    [10.0, 60.0, 110.0, 160.0, 210.0],
    [10.0, 43.3, 76.7, 110.0, 143.3],
    [10.0, 35.0, 60.0, 85.0, 110.0],
    [10.0, 30.0, 50.0, 70.0, 90.0],
    [130.0, 230.0, 330.0, 430.0, 530.0, 630.0],
    [80.0, 130.0, 180.0, 230.0, 280.0, 330.0],
    [63.3, 96.6, 130.0, 163.3, 196.6, 230.0],
    [55.0, 80.0, 105.0, 130.0, 155.0, 180.0],
    [50.0, 70.0, 90.0, 110.0, 130.0, 150.0]]

We set the initial weight of the synapse,

# set the initial weight of the synapse::
init_w = 0.5
syn_weights = []

and run the simulation for each pair of spike times (the weights are recorded for each simulation run):

for s_t_pre, s_t_post in zip(spike_times_pre, spike_times_post):
    nest.ResetKernel()
    nest.resolution = resolution

    # create one adaptive exponential integrate-and-fire neuron neuron with Clopath synapse:
    clopath_neuron = nest.Create("aeif_psc_delta_clopath", 1, clopath_neuron_params)

    # We need a parrot neuron for technical reasons since spike generators can only
    # be connected with static connections:
    parrot_neuron = nest.Create("parrot_neuron", 1)

    # create and connect spike generators:
    spike_gen_pre = nest.Create("spike_generator", {"spike_times": s_t_pre})
    nest.Connect(spike_gen_pre, parrot_neuron, syn_spec={"delay": resolution})

    spike_gen_post = nest.Create("spike_generator", {"spike_times": s_t_post})
    nest.Connect(spike_gen_post, clopath_neuron, syn_spec={"delay": resolution, "weight": 80.0})

    # create weight recorder:
    weightrecorder = nest.Create("weight_recorder")

    # create Clopath connection with weight recorder:
    nest.CopyModel("clopath_synapse", "clopath_synapse_rec", {"weight_recorder": weightrecorder})
    syn_dict = {"synapse_model": "clopath_synapse_rec", "weight": init_w, "delay": resolution}
    nest.Connect(parrot_neuron, clopath_neuron, syn_spec=syn_dict)

    # simulation:
    simulation_time = 10.0 + max(s_t_pre[-1], s_t_post[-1])
    nest.Simulate(simulation_time)

    # extract and save synaptic weights:
    weights = weightrecorder.get("events", "weights")
    syn_weights.append(weights[-1])

syn_weights = np.array(syn_weights)
# scaling of the weights so that they are comparable to Clopath et al (2010):
syn_weights = 100.0 * 15.0 * (syn_weights - init_w) / init_w + 100.0

Finally, we plot the results:

plt.figure(figsize=(4.5, 3.5))
plt.plot([10.0, 20.0, 30.0, 40.0, 50.0], syn_weights[5:], 
         lw=2.5, ls="-", label="pre-post pairing")
plt.plot([10.0, 20.0, 30.0, 40.0, 50.0], syn_weights[:5], 
         lw=2.5, ls="-", label="post-pre pairing")
plt.ylabel("normalized weight change")
plt.xlabel("firing rate of the presynaptic neuron [Hz]")
plt.legend()
plt.title(f"Synaptic weight using\nthe Clopath rule")
plt.tight_layout()
plt.show()

The resulting plot shows the synaptic weight change for different pairs of pre- and postsynaptic spike timings, plotted as a function of the firing rate of the presynaptic neuron. It shows how weight change depends not only on the order of spikes but also on the repetition frequency of the pairing protocol:

Synaptic weight change using the Clopath rule for different pairs of pre- and postsynaptic spike timings. The weight change is normalized to the initial weight and scaled to be comparable to Clopath et al. (2010)ꜛ. The two curves correspond to different temporal pairing orders and illustrate that the outcome depends jointly on spike timing, voltage dynamics, and repetition frequency.

The key observations in this plot are:

• pre-post pairing: As the firing rate of the presynaptic neuron increases, the synaptic weight change also increases. This indicates stronger synaptic potentiation at higher firing rates.
• post-pre pairing: The synaptic weight change increases similarly with higher firing rates but at lower rates compared to pre-post pairing. This shows a tendency for less potentiation or even depression in this pairing scenario for lower frequencies. However, for higher frequencies, the difference between pre-post and post-pre pairings diminishes.

The frequency of presynaptic spikes is crucial because it can significantly impact the degree of synaptic plasticity. Higher frequencies typically lead to more substantial changes in synaptic weights due to increased temporal overlap between pre- and postsynaptic spikes, which enhances the effects of spike-timing-dependent plasticity (STDP).

Network simulation of bidirectional connections

Next, we further investigate the Clopath rule by recapitulating the NEST tutorial “Clopath Rule: Large network simulation with bidirectional connections”ꜛ. In this simulation, we create a recurrent network of 10 excitatory and 3 inhibitory neurons connected with Clopath synapses, stimulated by Poisson generators. The goal is not to examine isolated spike pairs, but to see how voltage-based synaptic plasticity shapes connectivity within a small network of excitatory neurons embedded in an input driven circuit.

Let’s begin by importing the necessary libraries and setting the simulation parameters:

import random
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np
import nest

# set the verbosity of the NEST simulator:
nest.set_verbosity("M_WARNING")

# Set global properties for all plots
plt.rcParams.update({'font.size': 12})
plt.rcParams["axes.spines.top"]    = False
plt.rcParams["axes.spines.bottom"] = False
plt.rcParams["axes.spines.left"]   = False
plt.rcParams["axes.spines.right"]  = False


# set the simulation resolution and time:
simulation_time = 1.0e4
resolution = 0.1
delay = resolution
nest.ResetKernel()
nest.resolution = resolution

# for reproducibility:
np.random.seed(1)

Next, we define the parameters of the Clopath synapse and the neuron model and make connections between the neurons in the network. We create 10 excitatory and 3 inhibitory neurons and connect them with Clopath synapses. We also create Poisson generators to stimulate the network by 500 input neurons:

# poisson_generator parameters:
pg_A = 30.0     # amplitude of Gaussian
pg_sigma = 10.0 # std deviation

# create neurons and devices:
nrn_model = "aeif_psc_delta_clopath"
nrn_params = {
    "V_m": -30.6,           # [mV] membrane potential
    "g_L": 30.0,            # [nS] leak conductance
    "w": 0.0,               # [nS] adaptation conductance
    "tau_u_bar_plus": 7.0,  # [ms] time constant for u_bar_plus (LTP)
    "tau_u_bar_minus": 10.0,# [ms] time constant for u_bar_minus (LTD)
    "tau_w": 144.0,         # [ms] time constant for w
    "a": 4.0,               # [nS] subthreshold adaptation conductance
    "C_m": 281.0,           # [pF] membrane capacitance
    "Delta_T": 2.0,         # [mV] slope factor
    "V_peak": 20.0,         # [mV] spike cut-off
    "t_clamp": 2.0,         # [ms] clamping time
    "A_LTP": 8.0e-6,        # [nS] amplitude of LTP
    "A_LTD": 14.0e-6,       # [nS] amplitude of LTD
    "A_LTD_const": False,   
    "b": 0.0805,            # [nA] spike-triggered adaptation
    "u_ref_squared": 60.0**2# [nA^2] squared threshold for u
    }

pop_exc = nest.Create(nrn_model, 10, nrn_params) # create 10 excitatory neurons
pop_inh = nest.Create(nrn_model, 3, nrn_params)  # create 3 inhibitory neurons

pop_input = nest.Create("parrot_neuron", 500)    # helper neurons (for technical reasons)
pg = nest.Create("poisson_generator", 500)       # poisson generators; i.e., 500 input neurons
wr = nest.Create("weight_recorder")              # create a weight recorder

nest.Connect(pg, pop_input, "one_to_one", 
             {"synapse_model": "static_synapse", 
              "weight": 1.0, 
              "delay": delay})

nest.CopyModel("clopath_synapse", "clopath_input_to_exc", {"Wmax": 3.0})
conn_dict_input_to_exc = {"rule": "all_to_all"}
syn_dict_input_to_exc = {"synapse_model": "clopath_input_to_exc",
                         "weight": nest.random.uniform(0.5, 2.0),
                         "delay": delay}
nest.Connect(pop_input, pop_exc, conn_dict_input_to_exc, syn_dict_input_to_exc)

# create input->inh connections:
conn_dict_input_to_inh = {"rule": "all_to_all"}
syn_dict_input_to_inh = {"synapse_model": "static_synapse", "weight": nest.random.uniform(0.0, 0.5), "delay": delay}
nest.Connect(pop_input, pop_inh, conn_dict_input_to_inh, syn_dict_input_to_inh)

# create exc->exc connections:
nest.CopyModel("clopath_synapse", "clopath_exc_to_exc", {"Wmax": 0.75, "weight_recorder": wr})
syn_dict_exc_to_exc = {"synapse_model": "clopath_exc_to_exc", "weight": 0.25, "delay": delay}
conn_dict_exc_to_exc = {"rule": "all_to_all", "allow_autapses": False}
nest.Connect(pop_exc, pop_exc, conn_dict_exc_to_exc, syn_dict_exc_to_exc)

# create exc->inh connections:
syn_dict_exc_to_inh = {"synapse_model": "static_synapse", "weight": 1.0, "delay": delay}
conn_dict_exc_to_inh = {"rule": "fixed_indegree", "indegree": 8}
nest.Connect(pop_exc, pop_inh, conn_dict_exc_to_inh, syn_dict_exc_to_inh)

# create inh->exc connections:
syn_dict_inh_to_exc = {"synapse_model": "static_synapse", "weight": 1.0, "delay": delay}
conn_dict_inh_to_exc = {"rule": "fixed_outdegree", "outdegree": 6}
nest.Connect(pop_inh, pop_exc, conn_dict_inh_to_exc, syn_dict_inh_to_exc)

# set initial membrane potentials:
pop_exc.V_m = nest.random.normal(-60.0, 25.0)
pop_inh.V_m = nest.random.normal(-60.0, 25.0)

Finally, we simulate the network. We run the simulation in intervals of 100 ms and set the rates of the Poisson generators based on a Gaussian distribution:

# simulate the network:
sim_interval = 100.0
for i in range(int(simulation_time / sim_interval)):
    # set rates of poisson generators:
    rates = np.empty(500)
    # pg_mu will be randomly chosen out of 25,75,125,...,425,475
    pg_mu = 25 + random.randint(0, 9) * 50
    for j in range(500):
        rates[j] = pg_A * np.exp((-1 * (j - pg_mu) ** 2) / (2 * pg_sigma**2))
        pg[j].rate = rates[j] * 1.75
    nest.Simulate(sim_interval)

We then sort the synaptic weights according to the sender and reshape them into a matrix for visualization:

# sort weights according to sender and reshape:
exc_conns = nest.GetConnections(pop_exc, pop_exc)
exc_conns_senders = np.array(exc_conns.source)
exc_conns_targets = np.array(exc_conns.target)
exc_conns_weights = np.array(exc_conns.weight)
idx_array = np.argsort(exc_conns_senders)
targets = np.reshape(exc_conns_targets[idx_array], (10, 10 - 1))
weights = np.reshape(exc_conns_weights[idx_array], (10, 10 - 1))

# sort according to target:
for i, (trgs, ws) in enumerate(zip(targets, weights)):
    idx_array = np.argsort(trgs)
    weights[i] = ws[idx_array]

weight_matrix = np.zeros((10, 10))
tu9 = np.triu_indices_from(weights)
tl9 = np.tril_indices_from(weights, -1)
tu10 = np.triu_indices_from(weight_matrix, 1)
tl10 = np.tril_indices_from(weight_matrix, -1)
weight_matrix[tu10[0], tu10[1]] = weights[tu9[0], tu9[1]]
weight_matrix[tl10[0], tl10[1]] = weights[tl9[0], tl9[1]]

We also calculate the difference between the initial and final synaptic weights to assess the changes in synaptic weights over time:

# difference between initial and final value:
init_w_matrix = np.ones((10, 10)) * 0.25
init_w_matrix -= np.identity(10) * 0.25

Here are the plot commands to visualize the synaptic weight changes:

# plot synapse weights of the synapses within the excitatory population:
fig, ax = plt.subplots(figsize=(4.85, 4.5))
img = ax.imshow(weight_matrix - init_w_matrix, aspect='auto')

# create an axes on the right side of ax for the colorbar:
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(img, cax=cax)
cbar.ax.tick_params()
cbar.set_ticks(np.arange(-0.002, 0.0115, 0.002))

# adjust the labels and title positions:
ax.set_xlabel("to neuron")
ax.set_ylabel("from neuron")
ax.set_title("Change of synaptic weights\nbefore and after simulation")

# set x and y ticks:
xticklabels = ["1", "3", "5", "7", "9"]
ax.set_xticks([0, 2, 4, 6, 8])
ax.set_xticklabels(xticklabels)
yticklabels = ["1", "3", "5", "7", "9"]
ax.set_yticks([0, 2, 4, 6, 8])
ax.set_yticklabels(yticklabels)

plt.tight_layout()
plt.show()

Synaptic weight changes in a large network with 10 excitatory and 3 inhibitory receiving inputs from 500 Poisson generators. The plot shows the change in synaptic weights before and after the simulation. Rows indicate presynaptic neurons and columns indicate postsynaptic neurons. The synaptic weights are normalized to the initial weight.

The resulting matrix of weight changes shows the plasticity effects in the network, with positive values indicating potentiation and negative values indicating depression. The matrix is not diagonal-symmetric because the Clopath rule allows for bidirectional plasticity, meaning the changes in synaptic weights are not necessarily the same in both directions between any two neurons. The Clopath rule modifies synaptic weights based on the precise timing of spikes between the presynaptic and postsynaptic neurons. If neuron $i$ fires before neuron $j$, the weight change from $i$ to $j$ can be different from the weight change from $j$ to $i$, leading to an asymmetric weight matrix. And each synaptic connection updates its weight independently. The weight from neuron $i$ to neuron $j$ ($w_{i\rightarrow j}$) and from neuron $j$ to neuron $i$ ($w_{j\rightarrow i}$) are updated based on their own specific spike timing and neuronal activity. This independent update mechanism naturally leads to non-symmetric changes.

In biological neural networks, synaptic plasticity is inherently direction-dependent. The synaptic strength from one neuron to another can be modulated differently compared to the reverse direction, reflecting the biological processes of learning and memory encoding. The LTP and LTD processes depend on the relative timing of pre- and postsynaptic spikes. The temporal asymmetry in spike timing (e.g., pre-before-post vs. post-before-pre) results in different synaptic modifications for each direction.

Thus, the lack of diagonal symmetry in the weight change matrix is a direct result of the nature of the Clopath rule, which inherently supports bidirectional, independent, and asymmetric synaptic weight updates. This property allows for more complex and realistic modeling of synaptic plasticity, capturing the nuanced and direction-specific nature of biological synapses.

Conclusion

The Clopath rule is not just another STDP curve. It is a voltage-based plasticity model that explicitly links synaptic change to the postsynaptic depolarization state. This makes it far more suitable than simple pair-based STDP rules for modeling experiments in which firing rate, subthreshold voltage, and repeated spike pairing all influence the direction and magnitude of synaptic plasticity.

I think, its main strength lies in the combination of mechanistic interpretability and computational simplicity. The rule remains compact enough for large network simulations, yet rich enough to capture important biological features such as voltage dependence, frequency dependence, and homeostatic stabilization. For this reason, it has become a standard choice whenever synaptic plasticity is meant to depend on more than spike timing alone.

The complete code used in this blog post is available in this Github repositoryꜛ. The spike pairing and recurrent network examples can be adapted easily to explore how voltage-based plasticity shapes synaptic organization under different input regimes.

The complete code used in this blog post is available in this Github repositoryꜛ (clopath_spike_pairing.py and clopath_bidirectional_connections.py). Feel free to modify and expand upon it, and share your insights.

References

• C. Clopath, L. Büsing, E. Vasilaki, Wulfram Gerstner, Connectivity reflects coding: a model of voltage-based STDP with homeostasis, 2010, Nat Neurosci, 13(3), 344-352. doi: 10.1038/nn.2479ꜛ
• Claudia Clopath, Wulfram Gerstner, Voltage and spike timing interact in STDP - a unified model, 2010, Frontiers in Synaptic Neuroscience, Vol. n/a, Issue n/a, pages n/a, doi: 10.3389/fnsyn.2010.00025
• Voltage-based STDP synapse (Clopath et al. 2010) on ModelDBꜛ
• Wulfram Gerstner, Werner M. Kistler, Richard Naud, and Liam Paninski, Chapter 19 Synaptic Plasticity and Learning Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition, 2014, Cambridge University Press, ISBN: 978-1-107-06083-8, free online versionꜛ
• Jesper Sjöström, Wulfram Gerstner, Spike-timing dependent plasticity, 2010, Scholarpedia, 5(2):1362, doi: 10.4249/scholarpedia.1362ꜛ
• NEST’s tutorial “Clopath Rule: Spike pairing experiment”ꜛ
• NEST’s tutorial “Clopath Rule: Bidirectional connections”ꜛ
• NEST’s aeif_psc_delta_clopath model descriptionꜛ
• NEST’s clopath_synapse descriptionꜛ

Evolution of synaptic weights according to the Urbanczik-Senn plasticity model. We will discuss them in detail in the results part.

Core concepts

Unlike traditional STDP models, which rely solely on the relative timing of pre- and postsynaptic spikes, the Urbanczik-Senn plasticity model introduces a third factor: The local dendritic potential. The neuron model therefore consists of a somatic and a dendritic compartment. The role of the somatic compartment is to integrate inputs and generate spikes, while the dendritic compartment receives synaptic inputs. The local dendritic potential $V_d(t)$ predicts the somatic activity:

where $C_m$ is the membrane capacitance, $g_L$ is the somatic leak conductance, $E_L$ is the resting membrane potential, and $I_d(t)$ is the dendritic input current.

The somatic potential $V_s(t)$ is driven by the dendritic input $V_d$ and other synaptic inputs directly targeting the soma:

where $I_s(t)$ is the somatic input current, and $g_D$ is the dendritic coupling conductance. The model aims to minimize the discrepancy between the somatic firing rate $\phi(U)$ and the dendritic prediction $V_d^*$ of the somatic membrane potential:

The somatic firing rate $\phi(V_S)$ is given by:

where $\phi_{\text{max}}$ is the maximum rate, $k$ is the rate slope, $\beta$ is a parameter controlling the steepness, and $\theta$ is the threshold potential.

Synaptic weights are adjusted based on a local dendritic prediction error, which is the difference between the actual somatic spikes and the predicted firing rate from the dendritic potential:

Here, $\eta$ is the learning rate, $S(t)$ represents the somatic spike train, and $\frac{\partial V_d}{\partial w_i}$ is the contribution of synapse $i$ to the dendritic potential.

This plasticity rule is versatile and can support various learning paradigms:

• supervised learning: The somatic compartment receives a target signal that guides learning.
• unsupervised learning: The network generates its own teaching signals, promoting self-organization.
• reinforcement learning: The learning rate is modulated by a reward signal, enabling reinforcement learning.

The main advantage of this model is its ability to unify different learning paradigms under a single rule, driven by the dendritic prediction error. This rule is both biologically plausible and functionally powerful, offering a comprehensive framework for understanding and simulating synaptic plasticity in neural networks. It therefore enables more nuanced and robust learning dynamics compared to traditional STDP models:

• Predictive coding: The dendritic prediction error adjusts synaptic weights to better match somatic activity, embodying a form of predictive coding.
• Robust learning: By integrating dendritic potentials, the model captures more nuanced synaptic dynamics compared to traditional STDP, which only considers spike timings.
• Versatility: The model’s applicability to supervised, unsupervised, and reinforcement learning highlights its robustness and broad relevance to various neural processing tasks.

For further details, please refer to the original paper by Urbanczik and Senn (2014)ꜛ.

Simulation in NEST

In the following, we replicate the NEST tutorial “Weight adaptation according to the Urbanczik-Senn plasticity”ꜛ with some minor modifications. The simulation uses the pp_cond_exp_mc_urbanczik modelꜛ implemented in the NEST simulator. It consists of a two-compartment spiking point process neuron with conductance-based synapses and is capable to connect to a Urbanczik synapse urbanczik_synapseꜛ. The code reproduces the simulation results shown in figure 1B in Urbanczik’s and Senn’s original workꜛ and all simulation parameters are set accordingly except for the units: The simulation uses standard units instead of the unitless quantities used in the paper.

Let’s begin with importing the necessary libraries:

import os
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
import nest
import nest.raster_plot
# set the verbosity of the NEST simulator:
nest.set_verbosity("M_WARNING")

# Set global properties for all plots
plt.rcParams.update({'font.size': 12})
plt.rcParams["axes.spines.top"]    = False
plt.rcParams["axes.spines.bottom"] = False
plt.rcParams["axes.spines.left"]   = False
plt.rcParams["axes.spines.right"]  = False

Next, we need to define a couple of helper functions. First, we define two functions that generate inhibitory and excitatory inputs to the soma:

# define a function for the inhibitory input:
def g_inh(amplitude, t_start, t_end):
    """
    returns weights for the spike generator that drives the inhibitory
    somatic conductance.
    """
    return lambda t: np.piecewise(t, [(t >= t_start) & (t < t_end)], [amplitude, 0.0])

# define a function for the excitatory input:
def g_exc(amplitude, freq, offset, t_start, t_end):
    """
    returns weights for the spike generator that drives the excitatory
    somatic conductance.
    """
    return lambda t: np.piecewise(
        t, [(t >= t_start) & (t < t_end)], [lambda t: amplitude * np.sin(freq * t) + offset, 0.0]
    )

Then, we define a function that calculates the matching potential $U_M$ as a function of the somatic conductances $g_E$ and $g_I$:

This function computes the effective potential based on the contributions of excitatory and inhibitory synaptic inputs. It represents the equilibrium potential at the soma, considering the balance of excitatory and inhibitory inputs:

# define the matching potential:
def matching_potential(g_E, g_I, nrn_params):
    """
    returns the matching potential as a function of the somatic conductances.
    """
    E_E = nrn_params["soma"]["E_ex"]
    E_I = nrn_params["soma"]["E_in"]
    return (g_E * E_E + g_I * E_I) / (g_E + g_I)

Next, we define the dendritic prediction $V_d^*$ of the somatic membrane potential:

# define the dendritic prediction of the somatic membrane potential:
def V_w_star(V_w, nrn_params):
    """
    returns the dendritic prediction of the somatic membrane potential.
    """
    g_D = nrn_params["g_sp"]
    g_L = nrn_params["soma"]["g_L"]
    E_L = nrn_params["soma"]["E_L"]
    return (g_L * E_L + g_D * V_w) / (g_L + g_D)

We also need a function to calculate the rate function $\phi(U)$:

# define the rate function phi:
def phi(U, nrn_params):
    """
    rate function of the soma
    """
    phi_max = nrn_params["phi_max"]
    k = nrn_params["rate_slope"]
    beta = nrn_params["beta"]
    theta = nrn_params["theta"]
    return phi_max / (1.0 + k * np.exp(beta * (theta - U)))

and a function to calculate the derivative of the rate function $\phi$, i.e.,

# define the derivative of the rate function phi:
def h(U, nrn_params):
    """
    derivative of the rate function phi
    """
    k = nrn_params["rate_slope"]
    beta = nrn_params["beta"]
    theta = nrn_params["theta"]
    return 15.0 * beta / (1.0 + np.exp(-beta * (theta - U)) / k)

The derivative is needed to calculate the contribution of each synapse to the dendritic potential.

Now, we can start setting up the simulation. We first reset the NEST kernel and set the simulation parameters:

nest.ResetKernel()

# set simulation parameters:
n_pattern_rep     = 100  # number of repetitions of the spike pattern
pattern_duration  = 200.0
t_start           = 2.0 * pattern_duration
t_end             = n_pattern_rep * pattern_duration + t_start
simulation_time   = t_end + 2.0 * pattern_duration
n_rep_total       = int(np.around(simulation_time / pattern_duration))
resolution        = 0.1
nest.resolution   = resolution

Next, we set the neuron parameters, synapse parameters, and input parameters:

# set neuron parameters:
nrn_model = "pp_cond_exp_mc_urbanczik"
nrn_params = {
    "t_ref": 3.0,       # refractory period
    "g_sp": 600.0,      # soma-to-dendritic coupling conductance
    "soma": {
        "V_m": -70.0,   # initial value of V_m
        "C_m": 300.0,   # capacitance of membrane
        "E_L": -70.0,   # resting potential
        "g_L": 30.0,    # somatic leak conductance
        "E_ex": 0.0,    # resting potential for exc input
        "E_in": -75.0,  # resting potential for inh input
        "tau_syn_ex": 3.0,  # time constant of exc conductance
        "tau_syn_in": 3.0,  # time constant of inh conductance
    },
    "dendritic": {
        "V_m": -70.0,  # initial value of V_m
        "C_m": 300.0,  # capacitance of membrane
        "E_L": -70.0,  # resting potential
        "g_L": 30.0,   # dendritic leak conductance
        "tau_syn_ex": 3.0,  # time constant of exc input current
        "tau_syn_in": 3.0,  # time constant of inh input current
    },
    # set parameters of rate function:
    "phi_max": 0.15,    # max rate
    "rate_slope": 0.5,  # called 'k' in the paper
    "beta": 1.0 / 3.0,
    "theta": -55.0,
}

# set synapse params:
syns = nest.GetDefaults(nrn_model)["receptor_types"]
init_w = 0.3 * nrn_params["dendritic"]["C_m"]
syn_params = {
    "synapse_model": "urbanczik_synapse_wr",
    "receptor_type": syns["dendritic_exc"],
    "tau_Delta": 100.0,  # time constant of low pass filtering of the weight change
    "eta": 0.17,  # learning rate
    "weight": init_w,
    "Wmax": 4.5 * nrn_params["dendritic"]["C_m"],
    "delay": resolution,
}

# set somatic input:
ampl_exc = 0.016 * nrn_params["dendritic"]["C_m"] # amplitude of the excitatory input in nA
offset = 0.018 * nrn_params["dendritic"]["C_m"]   # offset of the excitatory input in nA
ampl_inh = 0.06 * nrn_params["dendritic"]["C_m"]  # amplitude of the inhibitory input in nA
freq = 2.0 / pattern_duration                     # frequency of the excitatory input in Hz
soma_exc_inp = g_exc(ampl_exc, 2.0 * np.pi * freq, offset, t_start, t_end) # excitatory input
soma_inh_inp = g_inh(ampl_inh, t_start, t_end)                              # inhibitory input

Then, we set the dendritic input by creating a spike pattern using Poisson generators. We record the spikes of a simulation of $n_{\text{pg}}$ Poisson generators and give the recorded spike times to spike generators:

# set dendritic input
n_pg = 200      # number of poisson generators
p_rate = 10.0  # rate in Hz
pgs = nest.Create("poisson_generator", n=n_pg, params={"rate": p_rate}) # poisson generators
prrt_nrns_pg = nest.Create("parrot_neuron", n_pg)                       # parrot neurons (for technical reasons)
nest.Connect(pgs, prrt_nrns_pg, {"rule": "one_to_one"})
spikerecorder = nest.Create("spike_recorder", n_pg)                # create the spike recorder
nest.Connect(prrt_nrns_pg, spikerecorder, {"rule": "one_to_one"})
nest.Simulate(pattern_duration)
t_srs = [ssr.get("events", "times") for ssr in spikerecorder]

After simulating the spike pattern, the spike times are stored in the variable t_srs and we need to reset the simulation kernel to start the actual simulation:

nest.ResetKernel()
nest.resolution = resolution

Now, we can create the neuron and the according spike generators:

# define neuron:
nrn = nest.Create(nrn_model, params=nrn_params) # create the Urbanczik neuron

# poisson generators are connected to parrot neurons which are connected to the mc neuron:
prrt_nrns = nest.Create("parrot_neuron", n_pg)

# create excitatory input to the soma:
spike_times_soma_inp = np.arange(resolution, simulation_time, resolution)
sg_soma_exc = nest.Create("spike_generator", 
                          params={"spike_times": spike_times_soma_inp, 
                                  "spike_weights": soma_exc_inp(spike_times_soma_inp)})
# create inhibitory input to the soma:
sg_soma_inh = nest.Create("spike_generator", 
                          params={"spike_times": spike_times_soma_inp, 
                                  "spike_weights": soma_inh_inp(spike_times_soma_inp)})

# create excitatory input to the dendrite:
sg_prox = nest.Create("spike_generator", n=n_pg)

We also create a multimeter for recording all parameters of the Urbanczik neuron, a weight recorder for recording the synaptic weights of the Urbanczik synapses, and another spike recorder for recording the spiking of the soma:

# create a multimeter for recording all parameters of the Urbanczik neuron:
rqs = nest.GetDefaults(nrn_model)["recordables"]
multimeter = nest.Create("multimeter", params={"record_from": rqs, "interval": 0.1})

# create a weight_recorder for recoding the synaptic weights of the Urbanczik synapses:
wr = nest.Create("weight_recorder")

# create another spike recorder for recording the spiking of the soma:
spikerecorder_soma = nest.Create("spike_recorder")

Finally, we connect all nodes:

# connect all nodes:
nest.Connect(sg_prox, prrt_nrns, {"rule": "one_to_one"})
nest.CopyModel("urbanczik_synapse", "urbanczik_synapse_wr", {"weight_recorder": wr[0]})
nest.Connect(prrt_nrns, nrn, syn_spec=syn_params)
nest.Connect(multimeter, nrn, syn_spec={"delay": 0.1})
nest.Connect(sg_soma_exc, nrn, 
             syn_spec={"receptor_type": syns["soma_exc"], 
                       "weight": 10.0 * resolution, 
                       "delay": resolution})
nest.Connect(sg_soma_inh, nrn, 
             syn_spec={"receptor_type": syns["soma_inh"], 
                       "weight": 10.0 * resolution, 
                       "delay": resolution})
nest.Connect(nrn, spikerecorder_soma)

and start the simulation, which is divided into intervals of the pattern duration:

# start the simulation, which is divided into intervals of the pattern duration:
for i in np.arange(n_rep_total):
    # Set the spike times of the pattern for each spike generator
    for sg, t_sp in zip(sg_prox, t_srs):
        nest.SetStatus(sg, {"spike_times": np.array(t_sp) + i * pattern_duration})
    nest.Simulate(pattern_duration)

After the simulation is completed, we read out the recorded data for plotting:

# read out devices for plotting:
# multimeter:
mm_events = multimeter.events
t = mm_events["times"]
V_s = mm_events["V_m.s"]
V_d = mm_events["V_m.p"]
V_d_star = V_w_star(V_d, nrn_params)
g_in = mm_events["g_in.s"]
g_ex = mm_events["g_ex.s"]
I_ex = mm_events["I_ex.p"]
I_in = mm_events["I_in.p"]
U_M = matching_potential(g_ex, g_in, nrn_params)

# weight recorder:
wr_events = wr.events
senders = wr_events["senders"]
targets = wr_events["targets"]
weights = wr_events["weights"]
times = wr_events["times"]

# spike recorder:
spike_times_soma = spikerecorder_soma.get("events", "times")

Here are the plot commands for plotting $V_s$ (membrane potential of the soma), $V_d$ (membrane potential of the dendrite), $V_d^*$ (dendritic prediction of the somatic membrane potential), $U_M$ (matching potential), somatic conductances $g_I$ and $g_E$, dendritic currents $I_{\text{in}}$ and $I_{\text{ex}}$, rates $\phi(U)$ and $\phi(V_d)$, and the rate derivative $h(V_d^*)$:

# plot the results:
lw = 1.0
fig1, (axA, axB, axC, axD, axE) = plt.subplots(5, 1, sharex=True, figsize=(8, 12))

# plot membrane potentials and matching potential:
axA.plot(t, V_s, lw=lw, label=r"$V_s$ (soma)", color="darkblue")
axA.plot(t, V_d, lw=lw, label=r"$V_d$ (dendrit)", color="deepskyblue")
axA.plot(t, V_d_star, lw=lw, label=r"$V_d^\ast$ (dendrit)", color="b", ls="--")
axA.plot(t, U_M, lw=lw, label=r"$U_M$ (soma)", color="r", ls="-", alpha=0.5)
axA.set_ylabel("membrane pot [mV]", )
axA.legend(loc="upper right")

# plot somatic conductances:
axB.plot(t, g_in, lw=lw, label=r"$g_I$", color="r")
axB.plot(t, g_ex, lw=lw, label=r"$g_E$", color="magenta")
axB.set_ylabel("somatic\nconductance [nS]")
axB.legend(loc="upper right")

# plot dendritic currents:
axC.plot(t, I_in, lw=lw, label=f"$I_$", color="r")
axC.plot(t, I_ex, lw=lw, label=f"$I_$", color="magenta")
axC.set_ylabel("dendritic\ncurrent [nA]")
axC.legend(loc="upper right")

# plot rates:
axD.plot(t, phi(V_s, nrn_params), lw=lw, label=r"$\phi(V_s)$", color="darkblue")
axD.plot(t, phi(V_d, nrn_params), lw=lw, label=r"$\phi(V_d)$", color="deepskyblue")
axD.plot(t, phi(V_d_star, nrn_params), lw=lw, label=r"$\phi(V_d^\ast)$", color="b", ls="--")
axD.plot(t, phi(V_s, nrn_params) - phi(V_d_star, nrn_params), lw=lw, 
         label=r"$\phi(V_s) - \phi(V_d^\ast)$", color="r", ls="-")
axD.plot(spike_times_soma, 0.15 * np.ones(len(spike_times_soma)), ".", 
         color="g", markersize=2, label="spike")
axD.legend(loc="upper right")

axE.plot(t, h(V_d_star, nrn_params), lw=lw, label=r"$h(V_d^\ast)$", color="g", ls="-")
axE.set_ylabel("rate derivative")
axE.legend(loc="upper right")
axE.set_xlim([0, 5000]) # we don't need to plot the whole simulation time

plt.tight_layout()
plt.show()

And here are the plot commands for plotting the evolution of synaptic weights:

# plot synaptic weights:
fig2, axA = plt.subplots(1, 1, figsize=(7, 4))
for i in np.arange(2, 200, 10):
    index = np.intersect1d(np.where(senders == i), np.where(targets == 1))
    if not len(index) == 0:
        axA.plot(times[index], weights[index], label="pg_{}".format(i - 2), lw=lw)

axA.set_title("Synaptic weights of Urbanczik synapses")
axA.set_xlabel("time [ms]")
axA.set_ylabel("weight")
axA.legend(fontsize=7, loc="upper right")
plt.tight_layout()
plt.show()

Results

In the following, we will take a look at the major results of the simulation, panel by panel. Note, that we limited the plots to the first 5000 ms (10000 ms for the synaptic weights) of the simulation to focus on the initial dynamics.

Membrane potentials and matching potential

First panel: membrane potentials of the soma ($V_S$, blue) and dendrite ($V_d$, cyan), the predicted dendritic potential ($V_d^*$, dashed blue) the matching potential ($U_M$, red).

The first panel shows the membrane potentials of the soma ($V_S$, blue) and dendrite ($V_d$, cyan), the predicted dendritic potential ($V_d^*$, dashed blue), and the matching potential ($U_M$, red). Let’s describe what we see here:

• Membrane potentials ($V_s$ and $V_d$): The somatic membrane potential ($V_s$) and dendritic membrane potential ($V_d$) show oscillatory behavior. This oscillation is driven by the excitatory and inhibitory inputs defined in your simulation parameters. The somatic membrane potential oscillates around a certain value, reflecting the integration of dendritic input and other somatic inputs.
• Predicted dendritic potential ($V_d^*$): This potential is calculated based on the somatic potential and the coupling conductance between soma and dendrite. It closely follows the dendritic potential ($V_d$), indicating that the model’s prediction aligns well with the actual dendritic activity.
• Matching potential ($U_M$): This potential is derived from the somatic excitatory and inhibitory conductances. It provides a reference for the neuron’s overall excitatory and inhibitory state.

Thus, the matching potential ($U_M$) aligns with the oscillations in the somatic and dendritic potentials, indicating that the neuron is balancing its excitatory and inhibitory inputs effectively. The predicted dendritic potential ($V_d^*$) closely tracking the actual dendritic potential ($V_d$) demonstrates the accuracy of the model’s prediction mechanism. This is crucial for the learning rule, as the synaptic adjustments depend on minimizing the discrepancy between the predicted and actual potentials. The somatic firing rate, which is not directly visible in this panel but is influenced by these potentials, will be regulated based on the predicted dendritic potential, ensuring that the neuron’s output remains consistent with its input patterns.

Somatic conductances and dendritic currents

Second panel: somatic conductances ($g_I$, red and $g_E$, magenta).

The second panel shows the somatic conductances ($g_I$, red, and $g_E$, magenta). $g_I$ starts at zero and quickly rises to a constant value, indicating a steady inhibitory input throughout the simulation. The excitatory conductance $g_E$ shows a sinusoidal pattern, oscillating in sync with the sinusoidal excitatory input applied. The steady inhibitory conductance ($g_I$) serves to balance the excitatory input and control the overall excitability of the neuron. The oscillatory excitatory conductance ($g_E$) directly reflects the applied sinusoidal input, showing that the synaptic inputs are effectively driving the somatic conductances as intended.

Dendritic currents

Third panel: dendritic currents ($I_{\text{in}}$, red, and $I_{\text{ex}}$, magenta).

The third panel shows the dendritic currents ($I_{\text{in}}$, red, and $I_{\text{ex}}$, magenta). The inhibitory current $I_{\text{in}}$ remains at zero, indicating that there is no inhibitory dendritic current being applied during the simulation. This is consistent with the settings used in the simulation, where only excitatory input was varied. The excitatory current $I_{\text{ex}}$ follows the sinusoidal pattern of the excitatory input, reflecting the synaptic input dynamics. The oscillatory excitatory current $I_{\text{ex}}$ oscillates with varying amplitude, showing periodic fluctuations over time. These fluctuations are influenced by the periodic excitatory input applied to the neuron, as per the settings of the simulation. This input drives the dendritic potential, which in turn affects the somatic potential and the overall activity of the neuron.

Overall, the oscillatory nature of the excitatory current reflects the sinusoidal excitatory input pattern defined in the simulation. The reduction in frequency over time may indicate a form of adaptation or a change in the neuron’s response to the input over time. The absence of inhibitory current suggests that the dynamics observed are primarily driven by excitatory inputs and their interaction with the neuron’s intrinsic properties.

Firing rates

Forth panel: fireing rates $\phi(U)$ (blue), $\phi(V_d)$ (light blue), $\phi(V_d^*)$ (blue dashed), $\phi(V_S) - \phi(V_d^*)$ (red), and the spikes (green dots).

The forth panel shows the different firing rates:

• $\phi(V_s)$ (blue): The somatic firing rate calculated from the somatic membrane potential $V_s$. It oscillates due to the synaptic inputs and the intrinsic dynamics of the neuron.
• $\phi(V_d)$ (light blue): The dendritic firing rate calculated from the dendritic membrane potential $V_d$. It follows a similar oscillatory pattern to the somatic firing rate but with some differences due to the separate dynamics of the dendritic compartment.
• $\phi(V_d^*)$ (blue dashed): The predicted dendritic firing rate, calculated from the predicted dendritic potential $V_d^*$. It also oscillates and aims to match the actual dendritic firing rate.
• $\phi(V_s) - \phi(V_d^*)$ (red): The difference between the somatic firing rate and the predicted dendritic firing rate. This discrepancy drives the synaptic plasticity according to the Urbanczik-Senn learning rule. Oscillations and deviations of this line from zero highlight periods where the model adjusts the synaptic weights to minimize this error.
• Spikes (green dots): The green dots along the top axis indicate the times at which the neuron fired action potentials (spikes). These spikes occur when the somatic membrane potential crosses a certain threshold, causing the neuron to emit an action potential.

Overall, this panel illustrates the dynamic interplay between the actual somatic and dendritic firing rates, the predicted dendritic firing rate, and the resulting discrepancies that drive synaptic adjustments. The goal of the learning rule is to minimize these discrepancies over time, leading to an adaptive and predictive neural response.

Rate derivative

Fifth panel: rate derivative $h(V_d^*)$.

The fifth panel shows the rate derivative $h(V_d^*)$ of the rate function $\phi(V_d^*)$. Initially, $h(V_d^*)$ starts at a high value around 5, indicating a steep rate of change in the firing probability. As time progresses, the value of $h(V_d^*)$ decreases, fluctuating between 1 and 4. This fluctuation shows the dynamic adjustment of the rate derivative in response to changes in the dendritic prediction potential. The overall trend shows an increase in fluctuations, indicating that the system is responding to varying synaptic inputs and adjusting the firing rate accordingly. The initial high value could be due to the neuron adapting rapidly to the initial conditions, after which it stabilizes and responds to the ongoing synaptic inputs

This rate derivative function plays a crucial role in the learning rule, as it affects the update of synaptic weights based on the dendritic prediction error. The fluctuations in $h(V_d^*)$ indicate active learning and adaptation processes in the neuron model, as it continuously adjusts to match the predicted somatic activity with the actual somatic firing rate.

Synaptic weights

Evolutions of synaptic weights of Urbanczik synapses.

This plot shows the evolution of synaptic weights for several synapses over time. Each colored line represents the weight of a synapse from a particular parrot neuron to the Urbanczik neuron. The weights exhibit different patterns of change, with some increasing significantly while others decrease or stabilize. The diversity in weight dynamics demonstrates the model’s capacity to differentiate between inputs based on their timing and interaction with the dendritic potential, leading to a self-organized pattern of synaptic strengths.

Overall interpretation

The plots collectively illustrate the dynamics of the Urbanczik-Senn plasticity model. The membrane potentials, conductances, and currents reflect the input-driven activity of the neuron. The firing rates and their discrepancies indicate the model’s predictive coding capabilities, where the dendritic compartment predicts somatic activity. The synaptic weights evolve based on the prediction errors, demonstrating the learning rule’s impact on synaptic plasticity.

Conclusion

The Urbanczik-Senn plasticity model offers a comprehensive framework for understanding synaptic plasticity in neural networks. By integrating dendritic prediction errors, the model unifies different learning paradigms under a single rule, enabling supervised, unsupervised, and reinforcement learning. The model’s predictive coding mechanism, robust learning dynamics, and versatility make it thus a powerful tool for simulating neural processing tasks and understanding the underlying mechanisms of synaptic plasticity.

The complete code used in this blog post is available in this Github repositoryꜛ (urbanczik_senn_plasticity.py). Feel free to modify and expand upon it, and share your insights.

References

• Robert Urbanczik, Walter Senn, Learning by the dendritic prediction of somatic spiking, 2014, Neuron, Vol. 81, Issue 3, pages 521-528, doi: 10.1016/j.neuron.2013.11.030ꜛ
• NEST’s tutorial “Weight adaptation according to the Urbanczik-Senn plasticity”ꜛ
• NEST’s pp_cond_exp_mc_urbanczik model descriptionꜛ
• NEST’s urbanczik_synapse synapse descriptionꜛ

In this post, we use nervosꜛ (Maskeen & Lashkare, 2025ꜛ) to implement a minimal two layer spiking neural network for pattern recognition on the MNIST dataset. We analyze how the network learns to classify digits through STDP and how the synaptic weights evolve during training. The figure shows the evolution of synaptic weights for a single output neuron over the course of training, illustrating how the network develops selectivity for certain input patterns corresponding to specific digit classes. We will further analyze the internal dynamics of the network, the learned receptive fields, and the classification performance on the test set. The code for this example is available in the Github repository mentioned at the end of this post.

I was curious and applied nervos to a classical pattern recognition task using MNIST digits. I basically replicated their original tutorial “MNIST Exampleꜛ” and then extended it with additional analyses. The goal was to understand in detail how a minimal two layer spiking network with purely local plasticity can self organize into digit selective neurons, how classification emerges, and what exactly is stored internally for later analysis. In this post, I summarize the mathematical formulation of the network, the learning rules, and the results I obtained.

Mathematics of nervos

Let’s first have a look at the mathematical formulation of the network architecture. In the following, we will walk through the main components of the model, including the network architecture, input encoding, neuron and synapse models, and the STDP learning rule. This is the best way to understand what the model does and how it learns, before we go ahead to the actual implementation and results.

Network architecture

The architecture implemented in nervos is deliberately minimal. It consists of a two layer spiking neural network without hidden layers:

• Input layer: we will use 784 neurons corresponding to the 28×28 pixels of an MNIST imageꜛ.
• Output layer: typically 60 or 80 neurons depending on the experiment (we will use 80 in our example).

Every input neuron is connected to every output neuron via an excitatory synapse. Thus, the weight matrix is

with entries $w_{ij}$ representing the synaptic strength from input neuron $j$ to output neuron $i$, with $N_{\text{out}}$ being the number of output neurons.

Competition in the output layer is implemented algorithmically. At each time step, the neuron with the highest membrane potential is identified. If this neuron crosses its adaptive threshold, it emits a spike and all other neurons are inhibited by resetting their potentials to an inhibitory level and placing them into a refractory state. There is no explicit lateral inhibitory connectivity matrix; instead, Winner Takes All dynamics are enforced by this global inhibition rule.

Samples from the MNIST datasetꜛ, which is used as input for the nervos SNN in our example below. Each image is 28x28 pixels, which corresponds to 784 input neurons in the network. The pixel intensities are converted to firing frequencies to generate spike trains for the input layer. The network learns to classify these images based on the spiking activity of the output layer neurons.

Input encoding

Each MNIST image is flattened into a vector $p \in [0,1]^{784}$. Pixel intensities are converted to firing frequencies according to

with typical values $f_{\max} = 70$ Hz and $f_{\min} = 5$ Hz.

For a presentation duration of $T$ discrete simulation steps, each input neuron generates a binary spike train

where $M_{ij} = 1$ if neuron $i$ fired at time step $j$.

Thus, the network receives a temporally structured, rate encoded spike pattern.

Neuron model: discrete integrate and fire dynamics with adaptive threshold

In practice, nervos does not implement a continuous time leaky integrate and fire differential equation with an explicit leak term. Instead, the neuron state is updated in discrete time steps. For each output neuron $i$ at time step $t$, the membrane potential is incremented by the weighted sum of incoming spikes

where $x_j(t)\in{0,1}$ is the presynaptic spike of input neuron $j$ at time step $t$.

After this synaptic integration step, nervos applies a discrete relaxation term whenever the neuron is above its resting potential $V_{\text{rest}}$:

with $\Delta_V=\texttt{spike_drop_rate}$. This term plays a stabilizing role similar to a leak, but it is not an exponential decay and it is not derived from a biophysical conductance model.

Each neuron also has a refractory mechanism implemented as a hard time step lock. After a neuron fires or is inhibited at time $t$, it is set to rest until

where $\tau_{\text{ref}}=\texttt{refractory_time}$. While $t < t_{\text{rest},i}$, the neuron does not integrate synaptic input.

The firing threshold is adaptive and implemented as an explicit state variable $\theta_i(t)$, initialized at a baseline value $\theta_0=\texttt{spike_threshold}$. Whenever a neuron fires, its adaptive threshold is increased additively

Additionally, when a neuron is above the resting potential, the adaptive threshold relaxes linearly toward the baseline by a fixed amount per time step

with $\Delta_\theta=\texttt{threshold_drop_rate}$. Thus, the adaptive threshold dynamics are discrete and piecewise linear rather than exponential with a time constant.

Finally, note that inhibition is implemented as a hard reset of the membrane potential to an inhibitory potential $V_{\text{inh}}=\texttt{inhibitory_potential}$ together with the same refractory lockout. This is a compressed, algorithmic representation of inhibition rather than an explicit inhibitory synapse conductance model.

Synapse models

nervos allows different synapse models, all normalized to

Three major models are implemented:

1. Ideal synapse: Continuous weight updates without quantization.
2. Linear finite state synapse: Uniform weight steps but with a finite number of discrete states.

3. Nonlinear memristor synaps: Based on experimental Pr$_{0.7}$Ca$_{0.3}$MnO$_3$ RRAM data. The $i$th weight state of an $n$ state synapse is modeled as

   \[\begin{aligned} w_i =& \quad w_{\max} - \frac{w_{\max} - w_{\min}}{1 - e^{-\nu}} \; \cdot \\ & \cdot \left[1 - \exp\left(-\nu\left(1 - \frac{i}{n}\right)\right)\right] \end{aligned}\]

   The parameter $\nu$ controls the curvature of the nonlinearity.

Initially, all synapses are set to $w = 1$ to facilitate early learning.

STDP learning

Weight changes are driven by an exponential STDP kernel $F(\Delta t)$ with

and

Note that nervos also supports alternative STDP kernels such as cosine, sinusoidal, and Gaussian depression-only variants. In the present analysis, however, we focus exclusively on the conventional exponential kernel defined above.

Synaptic updates are applied in a bounded, weight dependent manner. Let $w$ be the current weight and $w_{\min},w_{\max}$ the configured bounds. Define

Then the update is

with learning rate $\eta=\texttt{eta}$ and fixed exponent $\gamma=0.9$ in the current implementation. This makes potentiation and depression naturally saturate near the bounds.

STDP is applied only to synapses projecting onto a selected output neuron. In the default Winner Takes All mode (self.wta=True), synaptic updates are applied only for the winner neuron $k(t)$ at each time step where a spike event occurs. This is a strong form of competitive learning.

A further implementation detail is that if no presynaptic spike was observed in the configured past window for a given synapse at that postsynaptic event, nervos still applies a depression update with a randomly chosen negative $\Delta t$ value from a restricted range. This enforces ongoing weakening of synapses that are not consistently supported by correlated pre and post activity, thereby strengthening competition and sparsifying receptive fields.

Training and emergence of classification

During training, each image is presented as an input spike train for $T=\texttt{training_duration}$ discrete time steps. The network is simulated forward, and STDP updates are applied online whenever spiking events occur.

A key point is how nervos constructs the neuron to label association. In the current implementation, the neuron label map is updated online by directly assigning the true label of the current training sample to the neuron that was the maximally excited neuron when the last spike event occurred during that presentation. Denoting this neuron index by $k$, the update is

where $y$ is the true class label of the presented sample.

Thus, the label map is not computed via a separate majority vote over winners across the full training set. Instead, it is formed incrementally through repeated overwriting during training, and it stabilizes in practice because neurons that consistently win for a given class will repeatedly reassign themselves to that class.

At test time, classification proceeds without further plasticity. For a test spike train, the network computes winner event counts $c_i$ over the presentation window and returns

This readout is algorithmic and relies on the externally stored mapping from neurons to labels. It is therefore best interpreted as a minimal decision rule that extracts class predictions from the emergent winner selective dynamics of the trained output layer.

Synapse bounds and initialization

All synaptic weights are bounded by the configured limits $w_{\min}=\texttt{min_weight}$ and $w_{\max}=\texttt{max_weight}$. In the example configuration used below, these were set explicitly in the parameter dictionary, and the STDP update rule implements saturating weight dependence relative to these bounds.

In the current nervos implementation, the input to output synapse matrix is initialized as

that is all weights start at $w=1$. This choice accelerates early competition and receptive field formation, but it also means that early epochs can show very large weight norms and broad activation patterns before synaptic competition and bounded STDP drive weights toward more selective configurations.

What is stored internally

For detailed analysis, nervos stores (optionally):

• full weight matrix snapshots $W$ after each sample or epoch.
• spike raster matrices for each layer:
  $M^{(\text{layer})}_{\text{epoch}, \text{sample}}$
• learned neuron to label mapping
• synapse state trajectories for finite state models

From $W$, we can, e.g., compute receptive fields by reshaping

into $28 \times 28$ images, directly visualizing digit templates emerging in single neurons.

From spike rasters, we can compute:

• winner indices
• spike statistics
• L1 and L2 norms of synaptic vectors
• weight evolution curves

Thus, the framework allows a full dynamical and structural analysis of learning.

What nervos’ SNN does and does not achieve

While nervos is in my view a powerful tool to study STDP and synapse models, it is important to understand its strengths and limitations in the context of computational neuroscience and neuromorphic engineering.

Strengths

The main strengths of the nervos SNN are:

• fully local learning rule without global error backpropagation
• hardware aware synapse models including memristor nonlinearity
• very small architecture
• good performance under small training sets
• transparent internal dynamics

In small five class MNIST tasks, the authors report over 90% accuracy with conventional STDP and ideal synapses. When extended to all ten classes, accuracy drops to around 76%, which is still impressive for such a minimal architecture and purely local learning rule.

Limitations

The main limitations of the nervos SNN are:

• Two layer architecture only.
• Rate based input encoding, not true temporal coding.
• No deep hierarchical feature extraction.
• Classification relies on post hoc label assignment: It relies on an external label map that is built during training by repeatedly assigning the current sample label to the winning neuron, and is then used at test time to map the $\arg\max$ spike count neuron to a class label.
• Accuracy drops significantly when synapse states are strongly quantized.

Compared to fully biological plausible cortical microcircuits, the network is extremely simplified:

• no dendritic compartmentalization,
• no recurrent excitatory loops,
• no neuromodulation, and
• no reward signals.

Nevertheless, it provides a clean minimal platform to study STDP and hardware constraints.

Python example: Pattern Recognition on MNIST

Before we begin and for reproducibility, here is the environment setup that I have used for this example:

conda create -n nervos python=3.12 mamba -y
conda activate nervos
mamba install -y numpy matplotlib ipykernel requests
pip install nervos

I used nervos version 0.0.5, which is the latest version at the time of writing. The code is structured in a way that should be compatible with future versions, but some adjustments may be needed if the API changes significantly.

So, let’s start with the imports and some global settings for plotting:

import os
import numpy as np
import matplotlib.pyplot as plt
import nervos as nv

# set global properties for all plots:
plt.rcParams.update({'font.size': 12})
plt.rcParams["axes.spines.top"]    = False
plt.rcParams["axes.spines.bottom"] = False
plt.rcParams["axes.spines.left"]   = False
plt.rcParams["axes.spines.right"]  = False

Parameter setup

First, we define the parameters for the simulation. We will use 500 images of the MNIST training set for training and 150 images of the test set for testing. The training duration is set to 100 discrete time steps, which is sufficient for the network to process the input and generate spikes. We also set various parameters related to the neuron and synapse models, as well as the learning rates for STDP. We can also choose the number of classes we want to train on, e.g., 5 (MNIST subset) or 10 (full MNIST set). In this example, we will use 6 classes (digits 0 to 5) to keep the training time manageable while still demonstrating the learning capabilities of the network.

All parameters are stored in a Parameters object from the nervos library, which allows for easy access and modification throughout the code:

RESULTS_PATH = "figures"
os.makedirs(RESULTS_PATH, exist_ok=True)

# choose classes here:
CLASSES = list(range(6))   # choose any value between 1 and 10
identifier_name = f"{len(CLASSES)}classmnist"

p = nv.Parameters()

parameters_dict = {
    "training_images_amount": 500, # nervos is very memory hungry especially when m.get_spikeplots = True and m.get_weight_evolution = True (!); either use smaller numbers here, or set those to False below, or run it on a machine with high RAM
    "testing_images_amount": 150,
    "training_duration": 100, # discrete simulation time units
    "past_window": -10,
    "epochs": 3,  # note: after each epoch, the training set is not reshuffled, so the same images are presented in the same order.
    "image_size": [28, 28],
    "resting_potential": -70,
    "input_layer_size": 784,
    "output_layer_size": 80,
    "inhibitory_potential": -100,
    "spike_threshold": -55,
    "reset_potential": -90,
    "spike_drop_rate": 0.8,
    "threshold_drop_rate": 0.4,
    "min_weight": 1e-05,
    "max_weight": 1.0,
    "A_up": 0.8,
    "A_down": -0.3,
    "tau_up": 5,
    "tau_down": 5,
    "eta": 0.03,
    "min_frequency": 1,
    "max_frequency": 50,
    "refractory_time": 15,
    "tau_m": 10,
    "conductance": 10
}
for key, value in parameters_dict.items():
    setattr(p, key, value)

Note, that nervos is very memory hungry especially when m.get_spikeplots = True and m.get_weight_evolution = True (see below). Either use smaller numbers for the training and testing images (like ~500 and ~150, respectively), or set those settings to False, or run it on a machine with enough RAM.

Helper functions and class definitions

Next, we set up some functions which we need for the implementation of the MNIST_SNN class and for the analysis of the results.

We begin with the MNIST_SNN class which is a wrapper around the nervos SNN that handles data loading, training, and prediction. We also add a method to plot random samples from the dataset, which will be useful for visualizing the input spike trains before training the model. The plot_random_samples method takes care of aggregating the spike trains over time to reconstruct the original images for visualization, since the MNIST images are not stored as pixel values but only as spike trains in the model:

class MNIST_SNN(nv.Module):
    def __init__(self, parameters, identifier=None, classes=None, train_size=None, test_size=None, seed=None):
        super().__init__(parameters, identifier)

        # set default (5) if not provided
        if classes is None:
            classes = list(range(5))
        self.classes = list(classes)

        self.dataloader = nv.dataloader.MNISTLoader(parameters, classes=self.classes)

        # if you want the loader sizes to be controllable from outside:
        if train_size is None:
            train_size = getattr(parameters, "training_images_amount", 100)
        if test_size is None:
            test_size = getattr(parameters, "testing_images_amount", 20)

        self.X_train, self.Y_train = self.dataloader.dataloader(
            preprocess=True, pca=False, size=int(train_size), seed=seed)
        self.X_test, self.Y_test = self.dataloader.dataloader(
            preprocess=True, train=False, pca=False, size=int(test_size), seed=seed)

    def predict(self, un_processed_image, model_location):
        spike_train = np.array(self.dataloader.img2spiketrain(un_processed_image))
        synapses, neuron_label_map = self.load_model(model_location)
        return self.get_prediction(spike_train, synapses, neuron_label_map)

    def plot_random_samples(self, N=10, train=True, aggregate="sum", seed=None, cmap="hot_r", figsize=(10, 10)):
        """
        Plot N random MNIST samples from train or test set.

        Parameters
        ----------
        N : int
            Number of samples to plot.
        train : bool
            If True: use training set, else test set.
        aggregate : str
            "sum" or "mean" over time to reconstruct image from spike train.
        seed : int or None
            Optional seed for reproducibility.
            
            
        Note:
        -----
        nervos' dataloader directly returns spike trains for the MNIST images, 
        which we can visualize here before training the model. This also means,
        the MNIST images are not stored as pixel values in the model, but only 
        as spike trains. We therefore need to aggregate the spike trains over 
        time to reconstruct the original image for visualization.
        """

        if seed is not None:
            np.random.seed(seed)

        X = self.X_train if train else self.X_test
        Y = self.Y_train if train else self.Y_test

        indices = np.random.choice(len(X), size=min(N, len(X)), replace=False)

        cols = min(N, 5)
        rows = int(np.ceil(N / cols))

        plt.figure(figsize=figsize)

        for i, idx in enumerate(indices):
            spike_train = X[idx]  # shape: (784, T)

            if aggregate == "sum":
                img_vec = spike_train.sum(axis=1)
            elif aggregate == "mean":
                img_vec = spike_train.mean(axis=1)
            else:
                raise ValueError("aggregate must be 'sum' or 'mean'")

            img = img_vec.reshape(28, 28)

            plt.subplot(rows, cols, i + 1)
            plt.imshow(img, cmap=cmap, interpolation="nearest")
            plt.title(f"Label: {Y[idx]}")
            plt.axis("off")

        plt.tight_layout()

The next function is a direct replication of the original visualize_synapse function from the nervos tutorial, which visualizes the learned synaptic weights for each class. It aggregates the synaptic weights of all neurons that are assigned to the same class and reshapes them into 28x28 images to visualize the receptive fields learned by the network for each digit class:

def visualize_synapse(synapses, labels, cmap="hot_r", figsize=(10, 30), ncols=5):
    kk = 28
    labels = np.asarray(labels)

    classes = {i: np.zeros((kk, kk)) for i in np.unique(labels)}
    for idx in range(len(synapses)):
        classes[labels[idx]] += synapses[idx].reshape((kk, kk))

    class_keys = sorted(classes.keys())
    n_classes = len(class_keys)
    ncols = max(1, int(ncols))
    nrows = int(np.ceil(n_classes / ncols))

    plt.figure(figsize=figsize)
    for i, k in enumerate(class_keys, start=1):
        plt.subplot(nrows, ncols, i)
        plt.imshow(classes[k], cmap=cmap, interpolation="nearest")
        plt.title(f"{k}")
        plt.axis("off")

    plt.tight_layout()

To evaluate the performance of the model, we need to compute the accuracy on the test set. We use again a provided function from the nervos tutorial, called accuracy, which takes the trained model and the test classes as input, generates spike trains for the test images, and computes the predicted labels based on the output layer activity. It then compares the predicted labels to the true labels to calculate the overall accuracy of the model:

def accuracy(m2, classes, parameters_dict):
    loader = nv.dataloader.MNISTLoader(m2.parameters, classes=list(classes))
    spike_trains, labels = loader.dataloader(train=False, preprocess=True, 
      seed=123, size=parameters_dict["testing_images_amount"])

    t = 0
    c = 0
    preds = []
    print("Calculating Accuracy")
    for st, label in zip(spike_trains, labels):
        pred = m2.get_prediction(st)
        preds.append(pred)
        c += int(pred == label)
        t += 1
        print(f"\rTested {t} images", end="")
    print()
    print(c / t)
    return labels, preds

For a more detailed analysis of the model’s performance, we can compute the confusion matrix, which shows how many times each true class was predicted as each possible class. The following function confusion_matrix_np computes the confusion matrix which is used in the plot_confusion_matrix function to visualize it. The confusion matrix can be normalized to show percentages instead of raw counts, which can be helpful for interpretation:

def confusion_matrix_np(y_true, y_pred, labels):
    """
    Calculate the confusion matrix.
    Rows: true labels
    Cols: predicted labels
    """
    y_true = np.asarray(y_true, dtype=int)
    y_pred = np.asarray(y_pred, dtype=int)
    labels = np.asarray(labels, dtype=int)

    k = labels.size
    idx = {lab: i for i, lab in enumerate(labels)}
    C = np.zeros((k, k), dtype=int)

    for t, p in zip(y_true, y_pred):
        if (t in idx) and (p in idx):
            C[idx[t], idx[p]] += 1
    return C

def plot_confusion_matrix(C, labels, normalize=False, title=None, cmap="Greys"):
    """
    Plot confusion matrix. If normalize=True: row-normalize 
    (per true label).
    """
    C = np.asarray(C)
    if normalize:
        row_sums = C.sum(axis=1, keepdims=True).astype(float)
        row_sums[row_sums == 0.0] = 1.0
        M = C / row_sums
    else:
        M = C

    fig, ax = plt.subplots(figsize=(6, 5))
    im = ax.imshow(M, interpolation="nearest", cmap=cmap)
    cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)

    ax.set(
        xticks=np.arange(len(labels)),
        yticks=np.arange(len(labels)),
        xticklabels=[str(x) for x in labels],
        yticklabels=[str(x) for x in labels],
        xlabel="Predicted Labels",
        ylabel="True Labels",
        title=title if title is not None else ("Confusion matrix (normalized)" if normalize else "Confusion matrix"))

    # annotate:
    fmt = ".2f" if normalize else "d"
    thresh = M.max() * 0.6 if M.size else 0.0
    for i in range(M.shape[0]):
        for j in range(M.shape[1]):
            ax.text(
                j, i, format(M[i, j], fmt),
                ha="center", va="center",
                color="white" if M[i, j] > thresh else "black")

    plt.tight_layout()
    return fig, ax

The next function actually repeats the accuracy function but uses the confusion matrix to compute the overall accuracy. It takes the confusion matrix as input and calculates the accuracy as the trace (sum of diagonal elements) divided by the total number of samples. It also computes the recall for each class, which is the diagonal element divided by the sum of the corresponding row (true positives / total actual positives):

def accuracy_metrics(C):
    """
    From confusion matrix C (rows=true, cols=pred):
    accuracy.
    """
    C = np.asarray(C, dtype=float)
    total = C.sum()
    acc = np.trace(C) / total if total > 0 else np.nan

    row_sums = C.sum(axis=1)
    with np.errstate(divide="ignore", invalid="ignore"):
        recall = np.diag(C) / row_sums

    return acc, recall

An important aspect of analyzing the model’s performance is to look at the spike trains of the output neurons. The following rasterplot function creates a raster plot for the binary spike matrix, where each dot represents a spike from a neuron at a specific time step. The function also allows highlighting specific neurons (e.g., the epoch winner and the final winner) with different colors and sizes to visually distinguish them from the rest of the neurons:

def rasterplot(
    spike_train: np.ndarray,
    title: str = "raster",
    xlim=None,
    highlight_neuron_idx: int | None = None,
    highlight_color: str = "orange",
    highlight2_neuron_idx: int | None = None,
    highlight2_color: str = "magenta",
    base_color: str = "0.35",
    s_base: float = 2.0,
    s_highlight: float = 8.0,
    s_highlight2: float = 8.0):
    """
    Raster plot for a binary spike matrix with up to two highlighted neurons.

    spike_train: shape (N, T), entries 0/1
    highlight_neuron_idx: epoch wise winner (orange)
    highlight2_neuron_idx: final winner (magenta, retroactive)
    If both indices are equal, only one overlay is drawn (orange), but the legend
    label indicates "epoch winner = final winner".
    """
    spike_train = np.asarray(spike_train)
    if spike_train.ndim != 2:
        raise ValueError(f"spike_train must be 2D (N,T), got shape {spike_train.shape}")

    N, T = spike_train.shape
    ys, xs = np.where(spike_train == 1)

    plt.figure(figsize=(10, 4))
    plt.scatter(xs, ys, s=s_base, color=base_color, linewidths=0)

    handles = []

    def _overlay(idx: int, color: str, size: float, label: str):
        if not (0 <= idx < N):
            raise ValueError(f"highlight idx={idx} out of bounds for N={N}")
        xs_h = np.where(spike_train[idx] == 1)[0]
        if xs_h.size == 0:
            return None
        ys_h = np.full(xs_h.shape, idx, dtype=int)
        return plt.scatter(xs_h, ys_h, s=size, color=color, linewidths=0, label=label)

    j1 = int(highlight_neuron_idx) if highlight_neuron_idx is not None else None
    j2 = int(highlight2_neuron_idx) if highlight2_neuron_idx is not None else None
    same = (j1 is not None) and (j2 is not None) and (j1 == j2)

    # epoch winner (always, if provided)
    if j1 is not None:
        label1 = f"epoch winner = final winner idx {j1}" if same else f"epoch winner idx {j1}"
        h1 = _overlay(j1, highlight_color, s_highlight, label1)
        if h1 is not None:
            handles.append(h1)

    # final winner (only if provided and different)
    if (j2 is not None) and (not same):
        h2 = _overlay(j2, highlight2_color, s_highlight2, f"final winner idx {j2}")
        if h2 is not None:
            handles.append(h2)

    if handles:
        plt.legend(loc="best")

    if xlim is not None:
        plt.xlim(xlim)

    plt.xlabel("time step")
    plt.ylabel("neuron index")
    plt.grid(True, axis="x", linestyle="--", alpha=0.6)
    plt.title(title)
    plt.tight_layout()

To visualize the receptive fields of the output neurons, we define the following function plot_rf_of_neuron, which takes the synaptic weights of the first layer and a specific neuron index as input. It reshapes the weight vector of that neuron into a 28x28 image and visualizes it using a colormap. This allows us to see what kind of input pattern that neuron has learned to respond to:

def plot_rf_of_neuron(
    synapses_0: np.ndarray,
    neuron_idx: int,
    title: str = "",
    cmap: str = "viridis",
    figsize=(3.0, 3.0)) -> None:
    """
    Plot receptive field (weights reshaped to 28x28) of one output neuron.

    synapses_0: shape (n_out, 784)
    """
    w = np.asarray(synapses_0)[neuron_idx]  # (784,)
    img = w.reshape(28, 28)

    plt.figure(figsize=figsize)
    plt.imshow(img, cmap=cmap, interpolation="nearest")
    plt.title(title)
    plt.axis("off")
    plt.tight_layout()

The function plot_label_template visualizes a class specific weight template derived from the trained network. It selects all output neurons whose entry in neuron_label_map equals the given label and then averages their input weight vectors (or sums them, depending on mode). The resulting 784 dimensional vector is reshaped to 28×28 and plotted. This is therefore not an average MNIST image, but an average of learned synaptic weight patterns for that class. We use this function to visualize the emergent digit templates learned by the network for each class, which can be compared to individual receptive fields of single winner neurons as well as to the aggregated input spike representations of real MNIST samples, in order to assess how closely the learned synaptic structure aligns with the underlying data distribution:

def plot_label_template(
    synapses_0: np.ndarray,
    neuron_label_map: np.ndarray,
    label: int,
    title: str = "",
    cmap: str = "viridis",
    mode: str = "sum",  # "sum" or "mean"
    figsize=(3.0, 3.0)) -> None:
    """
    Plot label template: aggregate RFs of all neurons assigned to a given label.
    """
    synapses_0 = np.asarray(synapses_0)
    nlm = np.asarray(neuron_label_map)

    idx = np.where(nlm == int(label))[0]
    plt.figure(figsize=figsize)

    """ 
    idx.size indicates how many neurons are mapped to this label. 
    If idx.size == 0, it means no neuron is mapped to this label.
    """

    if idx.size == 0:
        plt.text(0.5, 0.5, f"No neurons mapped to label {label}", ha="center", va="center")
        plt.axis("off")
        plt.title(title)
        plt.tight_layout()
        return

    W = synapses_0[idx]  # (n_label_neurons, 784)
    if mode == "sum":
        img = W.sum(axis=0).reshape(28, 28)
    elif mode == "mean":
        img = W.mean(axis=0).reshape(28, 28)
    else:
        raise ValueError("mode must be 'sum' or 'mean'")

    plt.imshow(img, cmap=cmap, interpolation="nearest")
    plt.title(title + f" (neurons mapping: {idx.size}/{synapses_0.shape[0]})")
    plt.axis("off")
    plt.tight_layout()

The next functions are utility functions to extract the winner neuron index based on spike counts and to get the last weight snapshot for a specific sample and epoch. These functions will be useful for analyzing the evolution of the winner neuron’s receptive field over epochs:

def get_winner_neuron_idx(m, epoch, train_image_idx):
    spk_out = np.asarray(m.spikeplots[epoch][train_image_idx][-1])  # (n_out, T)
    spike_counts = spk_out.sum(axis=1)
    return int(np.argmax(spike_counts)), spike_counts

def get_last_weight_snapshot_for_sample(m, epoch, train_image_idx):
    """
    weight_evolution[epoch][sample] is a list of snapshots.
    Each snapshot has shape (n_out, n_in).
    """
    snapshots = m.weight_evolution[epoch][train_image_idx]
    if snapshots is None or len(snapshots) == 0:
        raise ValueError("No weight evolution stored. Set m.get_weight_evolution=True before training.")
    return np.asarray(snapshots[-1])  # (n_out, n_in)

def plot_winner_rf_evolution_over_epochs(
    m,
    train_image_idx,
    last_epoch=None,
    cmap="viridis",
    parameters_dict=None,
    nlm_final=None,):
    """
    For each epoch, pick the winner neuron by spike count
    (same definition as in your raster loop), then plot its RF from the
    epoch specific weight snapshot and compute summary metrics for that same winner.

    Notes:
      - The "winner" can change across epochs (that is the point).
      - If nlm_final is given, we also show map=... in the titles (final neuron label map).
    """
    if last_epoch is None:
        last_epoch = m.parameters.epochs - 1

    true_label = int(m.Y_train[train_image_idx])

    rfs = []
    norms_l1 = []
    norms_l2 = []
    means = []

    winner_idxs = []
    winner_counts = []
    winner_maps = []

    for ep in range(m.parameters.epochs):
        # epoch specific winner from spikes (same as your raster loop)
        spk_out = np.asarray(m.spikeplots[ep][train_image_idx][-1])  # (n_out, T)
        spike_counts = spk_out.sum(axis=1)
        winner_idx = int(np.argmax(spike_counts))
        winner_count = int(spike_counts[winner_idx])

        winner_map = None
        if nlm_final is not None and winner_idx < len(nlm_final):
            winner_map = int(nlm_final[winner_idx])

        # epoch specific weights (last snapshot for this sample at this epoch)
        W = get_last_weight_snapshot_for_sample(m, ep, train_image_idx)  # (n_out, 784)
        w = W[winner_idx]  # (784,)

        if w.size != 28 * 28:
            raise ValueError(
                f"Expected 784 weights for RF, got {w.size}. W shape is {W.shape}.")

        rfs.append(w.reshape(28, 28))
        norms_l1.append(np.sum(np.abs(w)))
        norms_l2.append(np.sqrt(np.sum(w**2)))
        means.append(np.mean(w))

        winner_idxs.append(winner_idx)
        winner_counts.append(winner_count)
        winner_maps.append(winner_map)

    # RF tiles
    n = len(rfs)
    fig, axes = plt.subplots(1, n, figsize=(3 * n, 3), squeeze=False)
    for ep in range(n):
        ax = axes[0, ep]
        ax.imshow(rfs[ep], cmap=cmap, interpolation="nearest")

        if winner_maps[ep] is None:
            ax.set_title(f"Epoch {ep}\nidx={winner_idxs[ep]}, spikes={winner_counts[ep]}")
        else:
            ax.set_title(
                f"Epoch {ep}\nidx={winner_idxs[ep]}, spikes={winner_counts[ep]}, map={winner_maps[ep]}"
            )

        ax.axis("off")

    plt.suptitle(
        f"Winner RF evolution for sample {train_image_idx}\ntrue={true_label}")
    plt.tight_layout()
    plt.savefig(os.path.join(RESULTS_PATH, f"winner_rf_evolution_sample{train_image_idx}_tiles.png"),dpi=200)
    plt.close()


    # summary metrics:
    fig, ax1 = plt.subplots(figsize=(6, 4))

    l1_line, = ax1.plot(norms_l1, marker="o", label="L1 norm")
    l2_line, = ax1.plot(norms_l2, marker="o", label="L2 norm")

    ax1.set_xlabel("Epoch")
    ax1.set_ylabel("L1 / L2 value")
    ax1.set_title(f"Weight summary for epoch wise winners (sample {train_image_idx})")
    # annotate eg at the L2 dots, which winner idx they correspond to:
    for ep in range(n):
        # ax1.annotate(f"winner\nidx {winner_idxs[ep]}", (ep, norms_l2[ep]), 
        #              textcoords="offset points", xytext=(0,10), ha='center', fontsize=8)
        ax1.annotate(f"winner\nidx: {winner_idxs[ep]}", (ep, np.max(norms_l1)*1.1), 
                     textcoords="offset points", xytext=(0,10), ha='center', fontsize=10,
                     va="top")
    ax1.set_ylim(bottom=0, top=np.max(norms_l1)*1.2)

    # secondary axis for mean weight:
    ax2 = ax1.twinx()
    mean_line, = ax2.plot(means, marker="o", linestyle="--", color="gray", label="Mean weight")
    ax2.set_ylabel("Mean weight", color="gray")
    ax2.tick_params(axis="y", labelcolor="gray")

    if parameters_dict is not None:
        ax2.set_ylim(bottom=parameters_dict["min_weight"], top=parameters_dict["max_weight"])
    else:
        ax2.set_ylim(bottom=0, top=1.01)

    # unified legend
    lines = [l1_line, l2_line, mean_line]
    labels = [line.get_label() for line in lines]
    ax1.legend(lines, labels, loc="best")
    plt.tight_layout()
    plt.savefig(os.path.join(RESULTS_PATH, f"winner_rf_evolution_sample{train_image_idx}_summary.png"), dpi=200)
    plt.close()