Design and Implementation of a Simulation Framework for a Bio–Neural Dust System

Abstract

This paper presents the development of a computer simulation framework, designed as a cost–effective and technically efficient alternative to experimental studies. The framework focuses on the Bio–Neural Dust System proposed in our previous works, which consists of two components: a light–emitting bio–nanosensor and an opsin–expressing genetically modified neuron. This innovative system holds significant potential for applications in neuroscience and biotechnology research. Programmed in Python, the framework provides researchers with a virtual tool to test and evaluate the Bio–Neural Dust System, enabling the prediction of outcomes for future in vivo experiments. This approach not only conserves resources, but also offers scientists a flexible and accessible means to investigate the complex interactions within the system prior to real–world applications. The framework’s adaptability and potential for diverse research applications highlight its importance in advancing the field of bio–nanotechnology.

1. Introduction

In the dynamic field of neuroscience, optogenetics stands as a pivotal research area, delving into the impact of light on genetically modified cells [1]. This discipline focuses primarily on genetically modified neurons that express opsins, proteins that are sensitive to light, thereby facilitating neuronal stimulation through specific wavelengths of light. Such research predominantly advances the fields of neuroscience and medical science. Traditional optogenetic methodologies entail the genetic modification of neurons to express opsins, which are then activated by light. This is commonly achieved by inserting a fiber optic filament into the brain or other regions rich in nerves, allowing for precise control over neuronal activity and, consequently, behavior, through light stimulation [2]. A significant example of this technique’s application is an experiment that induced predatory–like behavior in mice by stimulating their amygdala with blue light [3].

Despite its effectiveness, the invasive process of inserting fiber optic filaments raises considerable concerns, including the risks of infection and hemorrhage. From both ethical and practical standpoints, minimizing the invasiveness and associated risks of such brain–computer interface techniques is crucial. An innovative direction in this field is the elimination of these filaments in favor of a wireless approach to neuron stimulation. Recent advancements, such as Smith et al.’s (2024) development of an injectable fluorescent neural interface, have showcased non–invasive, cell–specific stimulation methods, further highlighting the potential for wireless optogenetic systems [4]. In 2018, a novel concept known as the Neural Dust System (NDS) was introduced [5]. The system comprised a microchip implanted onto a mouse’s sciatic nerve, which could be stimulated externally through ultrasonic waves [6]. While this represented a less invasive approach than previous methods, it still required surgical implantation.

This paper introduces a simulation framework for an advanced design of a Bio–Neural Dust System (BNDS), integrating a nanosphere transceiver that converts ultrasonic waves into light to stimulate opsin–expressing neurons. Notably, due to the size of the designed biological nanodevice, it can be suspended in a liquid medium (e.g., saline solution) and inserted into the brain cavity non–invasively, through injection, eliminating the need for inserting it surgically. Building upon our previous work in mathematical modeling, this study details the system’s components, simulation methodology, and results, along with future directions for development.

The rest of the paper is organized as follows. In Section 2, we delve into an extensive description of the Bio–Neural Dust System (BNDS), elucidating its innovative design and the intricate workings of its two primary components: the bio–optical transceiver and the opsin–expressing neuron. This section aims to provide readers with a thorough understanding of the BNDS architecture and its functionality. Section 3 is dedicated to outlining the methodology used in the development of the simulation framework. Here, we discuss the strategies and considerations involved in building a simulation model that accurately represents the BNDS. In Section 4, we shift our focus to the algorithms that were implemented within the simulation framework. Section 5 presents the outcomes of the simulations. In doing so, we aim to demonstrate the efficacy and accuracy of the simulation model in replicating the behavior of the BNDS components under various conditions. The section also ventures into future possibilities and improvements. Finally, Section 6 serves as the conclusion of the paper.

2. Bio–Neural Dust System: Overview

Figure 1 provides a representation of the system under investigation. As mentioned earlier, it comprises two primary components: the bio–optical transceiver and the opsin–expressing neuron. The bio–transceiver’s internal mechanisms will be discussed in detail later. However, at a fundamental level, the system operates through the following sequential steps:

• An external source emits an ultrasound at a predefined frequency and intensity.
• This ultrasonic wave penetrates the skin, skull, and brain tissue.
• Upon reception, the bio–transceiver processes this ultrasonic input, triggering its internal mechanisms to generate and emit blue light at a specific intensity (measured in mW/mm²).
• The opsins within the neuron membrane absorb this emitted blue light, leading to generating a current that potentially stimulates the neuron and may initiate an action potential.

Each of these steps and their underlying processes can be mathematically modeled, drawing on established principles in membrane physiology and neuroscience research. This section provides an overview of the high–level functions of each component and sub–component.

2.1. Bio–Transceiver

As delineated in reference [7], the bio–transceiver is conceptualized as a transparent nanosphere, composed of a hybrid bilayer and ranging in size from 1 μm to 1 mm (1000 μm). This nanosphere houses three key internal processes: the piezoelectric effect, calcium diffusion, and the bioluminescent reaction, each underpinned by distinct models but collectively functioning to convert ultrasonic waves into blue light.

2.1.1. Piezoelectric Effect

The piezoelectric effect is a phenomenon where certain materials convert mechanical force into electricity [8]. The M13 bacteriophage, a well–studied bacteriophage, exhibits this effect [9]. This component’s aim is to harness the M13 bacteriophage to absorb the ultrasonic wave and produce an electrical voltage for the subsequent component. This process involves anchoring the bacteriophage so that upon exposure to ultrasonic waves, it vibrates at its natural resonance frequency, generating a voltage [10]. This voltage is then transmitted through a nanowire to the next component.

2.1.2. Calcium Diffusion

This component involves a smooth endoplasmic reticulum (SER), which serves as a reservoir for Calcium ions (Ca²⁺). Typically found in muscle cells, the SER facilitates muscle contraction by releasing Ca²⁺ and reabsorbing it upon relaxation [11]. In this design, the SER is electrically connected to the M13 bacteriophage via a nanowire. When the M13 bacteriophage generates voltage, the SER responds by releasing Ca²⁺ into the nanosphere, where it interacts with the final component. The diffused ions exhibit Brownian motion along with Aequorin molecules within the sphere.

2.1.3. Bioluminescent Reaction

Aequorin, a photoprotein found in the jellyfish Aequorea Victoria, can bind with three Ca²⁺ ions to produce a photon of blue light (∼470 nm) [12]. In reference [7], Aequorin molecules are assumed to be spherical, about 100 ångstrom in size. As Ca²⁺ ions and Aequorin molecules move within the nanosphere, some ions bind to Aequorin, resulting in the emission of blue light.

The bio–transceiver’s processes described above form the foundation for the simulation framework detailed in the following sections. An understanding of these processes is crucial for the design and implementation of the simulation framework. The bioengineering feasibility of the key components is supported by insights from the existing literature. Different studies have addressed various aspects of the design. For instance, nanowire insulation can be achieved through established nanotechnology techniques, such as silicon dioxide or polymer–based coatings [13]. The stability of M13 bacteriophages, a critical component of our system, can be preserved by encapsulating them within polydimethylsiloxane (PDMS) [14].

The isolation of the sarcoplasmic reticulum (SR) from skeletal muscle is a well–documented process using differential centrifugation and sucrose density gradients, with evidence that isolated organelles can be sustained under appropriate conditions [15,16]. Additionally, the concentration of Aequorin within the system is feasible based on commercially available purified Aequorin products, which offer a solubility of up to 10 mg/mL. This corresponds to a maximum molar concentration of approximately 461 μM, calculated using the molecular weight of Aequorin (21.7 kDa), as follows:

$$\left[{\mathrm{Aq}}_{\max }\right]={\displaystyle \frac{\mathrm{Solubility} (\mathrm{g}/\mathrm{L})}{\mathrm{Molecular} \mathrm{Weight} (\mathrm{g}/\mathrm{mol})}}={\displaystyle \frac{10 \mathrm{g}/\mathrm{L}}{21700 \mathrm{g}/\mathrm{mol}}}\approx 0.000461 \mathrm{mol}/\mathrm{L}=461 \mu \mathrm{M}$$

(1)

Our simulations incorporate Aequorin concentrations within this feasible range to ensure practicality.

2.2. Opsin–Expressing Neuron

Neurons can be genetically altered to express opsins, light–sensitive proteins, in their cell membranes. When stimulated by light, these neurons can potentially fire an action potential. Optogenetics studies the induction of electrical activity in neurons using light.

2.2.1. Neurons and Membrane Potential

Neurons, the fundamental units of the brain and nervous system, are categorized into sensory, motor, and interneurons, each with distinct functions [17]. Sensory neurons process environmental stimuli, motor neurons control bodily movements, and interneurons facilitate neural signaling and contribute to higher–order processes such as cognition. The neuron structure includes a soma (cell body), an impermeable cell membrane, a signal–transmitting axon, and dendrites for receiving signals. The neuron’s electrical signaling is governed by its membrane potential, resulting from ion concentration differences across the cell membrane. Ion channels transport ions across this impermeable membrane, which is essential for electrical signaling. The membrane potential at rest is approximately –73 mV [18].

2.2.2. Opsins and Their Optogenetic Application

Opsins are light–sensitive proteins that function as ion transporters. The specific opsins in this system, Channelrhodopsin–2 (ChR2), are activated by blue light and transport Na⁺ ions into the cell, inducing a current that can trigger an action potential [19]. In optogenetics, neurons are genetically engineered to express ChR2, allowing light–induced stimulation of electrical activity. This application forms a cornerstone of optogenetic research, providing a novel method for manipulating neural behavior.

3. Framework Architecture

To create a software that effectively simulates the Bio–Neural system as outlined, it is essential to thoroughly understand and accurately model each component. An important aspect of this framework is to enable real–time visualization of the simulation, which is a key consideration in the implementation phase.

3.1. Bio–Transceiver

In order to make the simulations run in real–time, we need to dynamically track and calculate various variables, including the equivalent source voltage generated by the M13, SER membrane potential, concentrations of extracellular Ca²⁺ and Aequorin, and photons’ emission.

3.1.1. Piezoelectric Effect

Finite Element Analysis (FEA)–based studies in the literature [20,21,22] provided a dataset to understand the relationship between M13’s mechanical input and electrical voltage output at a macro scale. In our previous paper [10], we used FEA to study M13 nanogenerators, which means our work applies at the nanoscale. The method, discussed in detail in [10], captures the detailed physics of the nanogenerator, including its resonance frequency and voltage generation mechanisms.

In this study, we utilize a black box model derived from FEA to represent the M13 nanogenerator. While this approach simplifies the nanogenerator for system–level analysis, it remains grounded in the detailed physics obtained through FEA. Specifically, the input–output model used in this study reflects the behavior of an M13 nanogenerator in its cantilevered configuration at a fixed stimulation frequency ($f_{rc1}$ = 1.2985 MHz), matching the primary resonance frequency identified in the FEA. The model establishes a relationship between the intensity of ultrasonic waves and the resulting voltage output.

This dataset–driven approach offers flexibility for future studies. If the configuration of the nanogenerator changes, new FEA simulations would be required to update the black box model. To assist users in performing these updates, we have prepared a tutorial that outlines the process of obtaining the necessary parameters using FEA tools such as COMSOL Multiphysics. The tutorial includes step–by–step instructions on setting material properties, performing the analysis, and extracting the required input–output relationship. The tutorial has been made available in the project repository on GitHub. We should note that using open–source libraries to replicate the equivalent circuit proposed in [23,24,25,26] is another possible and more effective approach that we consider in our future work.

3.1.2. Calcium Diffusion

For simulating calcium diffusion, we focus on modeling the SER’s membrane potential fluctuation over time under external voltage. Various models exist for simulating membrane potential, each varying in complexity. In the current work, the Leaky Integrate and Fire (LIF) model [27] was chosen for its simplicity. This model iteratively calculates membrane potential over time increments ($d\tau$), using the following equations:

$$V_{steady-state}=V_{M13}+V_{rest}$$

(2)

$$V_t=V_{steady-state}+(V_{t-1}-V_{steady-state})e^{{\displaystyle \frac{-d\tau }{\tau }}}$$

(3)

where $V_t$ represents the current membrane potential, $V_{steady-state}$ is the steady–state voltage, $V_{M13}$ is the M13 bacteriophage’s source voltage, $V_{rest}$ is the membrane’s resting potential, $d\tau$ is the time step, and $\tau$ is the RC time constant of the membrane’s equivalent circuit. The membrane’s capacitance is constant, but its resistance varies with the number of calcium channels, as derived in [7].

To link the SER’s membrane potential to its extra/intracellular concentrations, necessary for calculating light emission, the Nernst equation is employed:

$$V={\displaystyle \frac{RT}{zF}}ln{\displaystyle \frac{\left[C{a}_{out}^{2+}\right]}{\left[C{a}_{in}^{2+}\right]}}$$

(4)

where R is the gas constant, T is the temperature in Kelvin, z is the valence of the Ca⁺² ion, F is the charge carried by one mole of electrons, $\left[C{a}_{out}^{2+}\right]$ is the Ca⁺² concentration outside the SER, $\left[C{a}_{in}^{2+}\right]$ is the Ca⁺² concentration inside the SER, and V is the SER’s membrane potential. This equation calculates the membrane’s reversal potential in relation to ion concentration gradients. As the membrane is permeable only to calcium ions, the resting potential equals the reversal potential. Applying an external voltage alters the membrane’s resting potential, which, in turn, changes the ion concentration gradient. This artificial alteration in resting potential induces a change in the concentration gradient. Using the LIF model, we can compute the new concentration gradient (extra/intra concentrations) based on the Nernst equation:

$$gradient={\displaystyle \frac{\left[C{a}_{out}^{2+}\right]}{\left[C{a}_{in}^{2+}\right]}}=e^{{\displaystyle \frac{V}{C}}}, with C={\displaystyle \frac{RT}{zF}}$$

(5)

Assuming a constant total calcium ion concentration within the nanosphere, we have:

$$\left[C{a}_{total}^{2+}\right]=\left[C{a}_{out}^{2+}\right]+\left[C{a}_{in}^{2+}\right]$$

(6)

This leads to two differential equations solved at each simulation interval:

$$\left[C{a}_{out}^{2+}\right]={\displaystyle \frac{\left[C{a}_{total}^{2+}\right]}{gradient+1}}$$

(7)

$$\left[C{a}_{in}^{2+}\right]=\left[C{a}_{total}^{2+}\right]-{\displaystyle \frac{\left[C{a}_{total}^{2+}\right]}{gradient+1}}$$

(8)

At the start of the simulation, $\left[C{a}_{total}^{2+}\right]$ is calculated and maintained as constant. Each time step includes measuring the membrane potential via the LIF model using Equation (2), determining the gradient, and then using differential equations to update the calcium ion concentrations.

3.1.3. Photon Emission

In our previous work [7], we presented our derived equations for estimating photon emission probability. The probability of a photon emitted from a single Aequorin molecule is:

$$P_{light}=1-P_{nolight}$$

(9)

$$P_{nolight}=(1+J_a)e^{-J_a\Delta T}$$

(10)

where $\Delta T$ is the time interval $d\tau$ and $J_a$ is the inward flux of Ca⁺² ions towards an Aequorin molecule.

$$J_a=4\pi D{r}_a\left[C{a}_{out}^{2+}\right]$$

(11)

where D is the Ca²⁺ diffusion coefficient and $r_a$ is the Aequorin’s radius. The rate of emitted photons per Aequorin is:

$$S_{LEM}={\displaystyle \frac{J_a}{3}}·P_{light}$$

(12)

The primary variable in these equations is $\left[C{a}_{out}^{2+}\right]$, updated by the Nernst equation model. As calcium ion concentrations change, Equation (11) calculates the rate of emitted photons per Aequorin. Multiplying this by the total number of Aequorin molecules gives the total emitted photon count. Finally, we calculate the equivalent blue light intensity:

$$lightintensity={\displaystyle \frac{\left[Ph\right]N_{A}hc}{\lambda }}$$

(13)

where $\left[Ph\right]$ is the total amount of emitted photons, $N_A$ is Avogadro’s number, h is Planck’s constant, c is the speed of light, and $\lambda$ is the wavelength of an emitted photon.

3.1.4. System Architecture Design

With the establishment of the required models, we are now equipped to execute the bio–transceiver simulation using Python.

Figure 2 illustrates the interconnectedness of each component within the system, demonstrating how they incorporate the defined models in their operations. Furthermore, the class diagram in Figure 3a outlines all the essential classes and scripts necessary for the operationalization of the simulation framework. The script designated as ’bubble_analysis.py’, positioned on the right in the diagram, is intended to be the primary script. It executes the algorithm by employing the mathematical models we have developed. On the other hand, ’bubble_3D.py’, located on the left, performs a similar function. However, it additionally offers a real–time 3D graphical representation of each model, enhancing the visualization aspect of the simulation. To further clarify the interactions between components in the system, we have added a Unified Modeling Language (UML) sequence diagram (Figure 3b) to complement the class diagram (Figure 3a). While the class diagram provides a static view of the relationships between components such as SmoothEndoRet and PhotonEmission, the sequence diagram illustrates the algorithmic flow and the relationships between variables and parameters used in the simulation framework. Specifically, it explains how data is processed within the framework, detailing the steps from ultrasonic input to updates in membrane potential, calcium ion dynamics, and photon emission, thereby highlighting the computational logic underpinning the system.

3.2. Opsin–Expressing Neuron

To develop a precise simulation of opsin–expressing neurons, two crucial elements were necessary. Firstly, a model was needed to accurately represent a neuron’s membrane potential, particularly its response to external stimuli. The second essential component was a method to simulate the photocurrent induced by the activation of Channelrhodopsin–2 proteins (ChR2), which serve as the stimuli triggering alterations in the membrane potential.

An extensive review of the relevant scientific literature revealed a range of mathematical models suitable for both the neuron’s membrane potential and the ChR2–induced photocurrent. These models vary in their levels of accuracy and computational efficiency. A comparative analysis of the models for each component was conducted to identify those that offer an optimal balance between precision and computational resource demands. The Izhikevich model was initially selected for its balance between computational efficiency and accuracy [28], making it well–suited for simulating the neuron’s membrane potential in early stages of the framework’s development. However, as our simulations progressed, it became evident that greater precision was required to accurately capture the photocurrent induced by Channelrhodopsin–2 proteins. Consequently, the Four–State Model was adopted for its ability to provide a more detailed representation of the ChR2 photocurrent dynamics while maintaining reasonable computational efficiency [29].

System Architecture Design

The development of the simulation framework for opsin–expressing neurons required the integration of two selected models, while ensuring the design remained adaptable for potential future replacements or updates to these models. In the proposed system configuration, depicted in Figure 4, the process starts with the channelrhodopsin component, which receives light intensity (irradiance) as input. This component then utilizes the Four–State Model [29] to calculate the resultant photocurrent. This generated current is subsequently relayed to the neuron component. Here, the Izhikevich model [28] is applied to assess the membrane potential, determining the occurrence of an action potential.

The system’s architecture design is outlined in a simplified class diagram in Figure 5a and clarified further with the sequence diagram in Figure 5b. Central to this design is the neuron class, which encompasses the overall system functionality by holding instances of both the Izhikevich model and the Channelrhodopsin–2 implementation. The nextTimeStep is responsible for calculating the membrane potential and photocurrent at subsequent time intervals, drawing on the functionalities of these two model instances. The Channelrhodopsin class, designed for further abstraction, is initialized with parameters consistent across various ChR2 models. It houses an instance of the FourStateModel class, which is called upon in its nextTimeStep for photocurrent computations. The FourStateModel class is set up with crucial parameters $O_1,O_2,C_1,C_2,$ and $p,$ using initial values detailed in the github repository. Additional necessary parameters are also established or computed within the getPhotocurrent function. This function is tasked with implementing the Four State Model equations to calculate the photocurrent of ChR2 molecules, taking into account variables such as gChR2 and EChR2, for a given irradiance, wavelength, and holding potential (V) over the next time interval. Conversely, the Izhikevich class initializes with v representing the resting potential, u as the product of b and v, and key parameters a, b, c, and d. The updateMembranePotential in this class executes the Izhikevich model equations, computing the neuron’s voltage at a specific time, given the input current.

This architectural design offers the necessary abstraction for the model components, ensuring that the system’s functionality is not heavily dependent on the specific computational mechanisms of the models. It allows for the replacement of the Izhikevich class with any other model capable of calculating membrane potential and contains updateMembranePotential $(time,current)$. Similarly, the FourStateModel class can be substituted with any alternative ChR2 model that includes a getPhotocurrent $(irradiance, V, dt, wavelength, g\_ChR2$, $E\_ChR2)$, facilitating future modifications or enhancements to the system.

4. Implementation

The simulation framework, along with its supporting models, has been developed using Python. The complete source code is available in the GitHub repository at: https://github.com/Arash-Azarnoush/BioNeuralDust-PublicRelease, accessed on 06 December 2023. Within this repository, distinct folders contain the scripts for each of the two primary components.

Beyond just implementing the mathematical models detailed previously, this framework also enables the visualization of a 3D representation of these models. Using Python for development, VPython 7, an open–source module for 3D graphics, was utilized for this purpose.

4.1. Bio–Transceiver

The /nano–bubble/directory in the repository contains all scripts related to the simulation framework. A brief overview of each script and the main script’s algorithm is as follows:

4.1.1. Main Scripts: bubble_analysis.py and bubble_3D.py

-bubble_analysis.py: This script aggregates and executes all the necessary classes and functions to run the simulation framework’s mathematical models.

-bubble_3D.py: Similar to bubble_analysis.py, this script also facilitates the real-time 3D rendering of the simulation.

4.1.2. Piezoelectric Effect

M13.py script includes the M13Phage (fileName) class, used to emulate the piezoelectric effect of the M13 bacteriophage. The fileName parameter specifies the dataset path from the finite element field analysis, available in the repository as M13_Voltage_v3.xlsx.

-setUltrasound (intensity, channels): This function takes the input intensity from the user and updates the source voltage using a dataset lookup.

-getVoltage (): Returns the current source voltage emitted by the bacteriophage. M13_3D.py inherits from the M13 class, adding VPython 3D geometry for visualization. It is used by bubble_3D.py.

4.1.3. Calcium Diffusion

ser.py contains the SmoothEndoRet () class, which abstracts the SER membrane and its changing resting potential using the LIF model. It tracks variables like $V_t$, $V_{steady-state}$, $V_{M13}$, $V_{rest}$, and the RC time constant $\tau$.

-updateVoltage VM13 updates the new source voltage VM13 and VSteadysSate based on Equation (1).

-nextTimeStep () computes the membrane potential Vt for the next time step using Equation (2). calcium.py features the CalciumCluster () class, implementing the Nernst equation to calculate intra– and extracellular Ca²⁺ concentrations. -getPoutMax () determines the maximum possible extracellular Ca²⁺ concentration.

-updateConcentration (Vt) uses Equations (4)–(7) to update Ca²⁺ concentrations, taking $V_t$ from ser.nextTimeStep () as input.

4.1.4. Photon Emission

photon_emission.py includes the PhotonEmission () class, which calculates the total blue light intensity from emitted photons based on the extracellular Ca²⁺ concentration.

-getLightIntensity (concentration, dt) accepts an extracellular Ca²⁺ concentration and computes Equations (8) to (12) to yield the resultant blue light intensity.

This comprehensive framework ensures a modular and flexible approach to simulating the Bio–Neural Dust system, with each component accurately modeled and visually represented.

4.1.5. Algorithm

The primary scripts of our Python–based simulation framework follow the Algorithm 1 outlined below:

Algorithm 1 Simulation of the bio–transceiver

Initialize the calcium model with CalciumCluster(). Define function zap(input_soundwave, channels): Set ultrasound parameters in bacteriophage with setUltrasound(). Retrieve voltage from bacteriophage with getVoltage(). Instantiate SmoothEndoRet, resting, and max voltages. Update SER voltage with updateVoltage(). Compute getLightIntensity(). Initialize Vt_array, Pt_out_array, light_array. Set initial values: Vt_array, Pt_out_array, light_array. For each time step t in S_array:

Update membrane potential with ser’s nextTimeStep().

Update updateConcentrations().

Compute getLightIntensity().

Store computed Vt_array, Pt_out_array, and light_array.

Return tuple of Vt_array, Pt_out_array, and light_array.

4.2. Opsin–Expressing Neuron

The scripts that are crucial for simulating genetically modified neurons are in the /neuron_sim/ directory of the project repository. This section aims to provide a comprehensive overview of the implementation process and the outcomes of the simulations.

4.2.1. System Logic

The computations of the neuron’s changing membrane potential are managed by the Izhikevich class found in /neuron_sim/izhikevich.py. This class implements the model’s equations, previously discussed, within its updateMembranePotential. Referring to the pseudocode:

1: function updateMembranePotential(t, I) 2:     $dv\leftarrow {(0.4·v)}^2+5·v+140-u+I$ 3:     $du\leftarrow a·(b·v-u)$ 4:     $v\leftarrow v+(dv·dt)$ 5:     $u\leftarrow u+(du·dt)$ 6:     if$v > 30$then 7:         $v\leftarrow c$ 8:         $u\leftarrow u+d$ 9:     end if 10: end function

• Line 2 calculates the change in voltage (v) due to an input current (I).
• Line 3 calculates the change in the recovery variable (u).
• Lines 5 and 6 update v and u for the time interval (dt).
• Lines 8–10 handle the variables’ reset post–spike.

The FourStateModel class in /neuron_sim/four_state_model.py includes the getPhotocurrent, which is crucial for determining the light–induced current of opsins:

1: function getPhotocurrent(i, V, $dt$, $wavelength$, $g_{ChR2}$, $E_{ChR2}$) 2:     $So\leftarrow 0.5·(1+tanh(120·(100·i-0.1)\left)\right)$ 3:     $dP\leftarrow (So-p)/T_{ChR2}$ 4:     $dO1\leftarrow ({ϵ}_{p1}·F·p·c1)+({ϵ}_{21}·o2)-(G_{d1}·o1)-({ϵ}_{12}·o1)$ 5:     $dO2\leftarrow ({ϵ}_{p2}·F·p·c2)+({ϵ}_{12}·o1)-(G_{d2}·o2)-({ϵ}_{21}·o2)$ 6:     $dC1\leftarrow (G_r·c2)+(G_{d1}·o1)-({ϵ}_{p1}·F·p·c1)$ 7:     $dC2\leftarrow (G_{d2}·o2)-({ϵ}_{p2}·F·p·c2)-(G_r·c2)$ 8:     $p\leftarrow p+dP·dt$ 9:     $c1\leftarrow c1+dC1·dt$ 10:     $o1\leftarrow o1+dO1·dt$ 11:     $o2\leftarrow o2+dO2·dt$ 12:     $c2\leftarrow c2+dC2·dt$ 13:     $I_{ChR2}\leftarrow g_{ChR2}·((10.6408-14.6408·exp(-V/42.7671))/V·(o1+\gamma ·o2)·(V-E_{ChR2})$ 14:     return  $I_{ChR2}$ 15: end function

• It implements equations for updating key parameters ($p,C_1,O_1,O_2,C_2$) over time (dt).
• It computes the induced current ($I_{ChR2}$) based on these parameters.

The Channelrhodopsin class in /neuron_sim/channelrhodopsin.py primarily acts as an abstraction layer. Its nextTimeStep invokes the FourStateModel instance for photocurrent computation. Similarly, the Neuron class in /neuron_sim/neuron.py is also an abstraction layer, combining membrane potential and photocurrent components. Its nextTimeStep calls on instances of Channelrhodopsin and Izhikevich for computations.

1: function nextTimeStep(t, $irradiance$) 2:     $opsinModel.nextTimeStep($ 3:         $membraneModel.getMembranePotential\left(\right),$ 4:         $irradiance)$ 5:     $membraneModel.updateMembranePotential($ 6:         $t,opsinModel.getI\left(\right))$ 7:     $apVisual(membraneModel.isSpiking(\left)\right)$ 8:     $opsinStateVisual\left(\right)$ 9: end function

The core functionality of the simulation, including the initialization and execution of various components, is governed by the main program located in /neuron_sim/main.py. This program ensures that the simulation runs continuously until it is manually terminated. As illustrated in the loop, a perpetual while true loop (line 1) maintains the simulation’s active status, even after the designated runtime (sim_time) concludes. This design facilitates the restarting of the simulation without having to restart the entire program. Within this loop, if the current simulation time (t) has not yet reached its scheduled endpoint (sim_time) and the simulation has not been paused (running is true), the program iterates through lines 4–17. At each iteration, the neuron instance is invoked (line 4) to calculate the membrane potential and photocurrent. This loop concludes with an incrementation of the simulation time (t) by the predefined time step (dt) on line 17, thus progressing the simulation.

1: while True do 2:     while  $t<sim\_time$  do 3:         if running then 4:             $neuron.nextTimeStep(t,getIrradiance(\left)\right)$ 5:             if  $getIrradiance\left(\right)>0$  then 6:                 $light\_source.visible\leftarrow True$ 7:             else 8:                 $light\_source.visible\leftarrow False$ 9:             end if 10:             $update(V\_display.text)$ 11:             $ChR2\_plot.plot\left(\right[t,neuron.getOpsin$ 12:                 $Model\left(\right).getI\left(\right)\left]\right)$ 13:             $neuron\_plot.plot\left(\right[t,neuron.getMembrane$ 14:                 $Model\left(\right).getMembranePotential\left(\right)\left]\right)$ 15:             $light\_plot.plot\left(\right[t,getIrradiance\left(\right)\left]\right)$ 16:             $t\leftarrow t+dt$ 17:         end if 18:     end while 19: end while

4.2.2. 3D Neuron Model

A significant effort went into creating a visually accurate 3D model of a neuron within the Neuron class:

• The soma is represented as a simple sphere.
• The axon, including its myelin sheath, is depicted using a series of ellipsoids forming a curved path.
• The axon terminals are modeled as curved branches with cone–like ends.
• The dendrites are the most complex, designed with selective randomness to mimic their natural structure.

4.2.3. Visual Effects

The simulation incorporates various visual effects to signify different events:

• A blue light illuminates the neuron model during optical stimulation.
• A brief flash at the axon terminals indicates the firing of an action potential.
  The proportion of activated ChR2 molecules is visually represented by coloring a fraction of the opsins blue (activated state) and the rest red (closed state).

4.2.4. User Interface

The simulation’s user interface is designed for a seamless user experience, featuring the following:

• Start/stop and reset buttons for controlling the simulation;
• A real–time display of simulation time and membrane potential;
• Graphs for photocurrent, membrane potential, and light intensity, updated in real–time;
• A control panel for adjusting simulation parameters, with changes immediately affecting the simulation.

The combination of these elements results in a highly functional and user–friendly simulation interface (Figure 6).

5. Results

5.1. Bio–Transceiver

The execution of bubble_analysis.py yields several graphs, conducted under specific initial conditions: an extracellular concentration of 0.2 millimoles, a resting potential $V_{rest}$ of –50 mV, a maximum potential $V_{max}$ of –10 mV, and a diffusion coefficient of 10⁻⁸. When operating the simulation with its 3D variant, bubble_3D.py, the program is launched in a browser window. This interface allows the user to adjust the ultrasonic intensity using a slider positioned at the top right of the window. Initiating the script enables the user to observe the changes in membrane potential, extracellular concentration, and blue light intensity in real–time, as showcased in Figure 6. The simulation visually represents the emission of blue light from the Aequorin as the concentration of calcium ions increases. Additionally, the 3D simulation offers interactive features, such as the capability to navigate around the canvas, enhancing the user experience. While bubble_3D.py performs identical computations to bubble_analysis.py, it tends to operate at a slower pace. This reduction in speed is attributed to the computational demands of rendering numerous spheres and objects during the simulation.

Overall, the results of the simulation are in line with the anticipated outcomes. However, there is a noticeable point of consideration regarding the initial extracellular concentration required to achieve the light intensities depicted in Figure 7. This discrepancy may stem from the underestimation of the lower limit of light emission intensity. Further scrutiny into the simulation parameters and the refinement of the light emission modeling could help resolve this issue.

5.2. Opsin–Expressing Neuron

Figure 8 shows the designed simulation’s user interface with the neuron’s 3D model and ChR2 molecules. The real–time display of the membrane potential shows that the neuron is at rest and therefore, ChR2 molecules are coloured in red (closed state). The figure also shows the control panel for adjusting simulation parameters. The graphs for photocurrent, membrane potential, and light intensity generated by the interface are shown in Figure 9.

The development of our system involved a meticulous validation process, aiming to ensure the accuracy of its results. This process entailed a comparison of the simulation’s outputs with established experimental data, followed by iterative adjustments to the model until it yielded satisfactory outcomes. Our primary reference for validation was the data from [29], which significantly informed our model’s design.

Initially, the model was validated against experimental data involving the photocurrent induced by ChR2 in response to a 500 ms light pulse at a 470 nm wavelength. This comparison was made across three different combinations of holding potential and irradiance (–20 mV & 5.5 mW/mm², –80 mV & 0.34 mW/mm², –80 mV & 5.5 mW/mm²). Figure 10 shows a comparison between (a) the experimental results in [29] and (b) our simulation outputs. The red line in Figure 10 and Figure 11 represents our simulation results overlaid transparently on the experimental data to visually demonstrate the high accuracy of the simulation in replicating the ChR2 photocurrent dynamics, despite minor discrepancies. This comparison demonstrates the high accuracy of the simulation framework in replicating experimental trends for ChR2 photocurrents. The small fluctuations observed can be attributed to differences in experimental setups and parameterization rather than limitations of the model. Experimental data naturally include variability due to stochastic biological processes, while the simulation operates under controlled conditions, emphasizing core dynamics. Such variabilities contributing to these differences include slight discrepancies in light absorption and scattering, actual neuron–specific opsin expression levels, and dynamic photocurrent responses. These minimal discrepancies highlight the strength of the framework in capturing essential dynamics while leaving room for refinement in experimental alignment. The simulation’s fidelity was further corroborated by comparing its output with the results of the same model used in the reference study [29], as shown in Figure 11. This validation involves a set of experimental data related to the photocurrent response to two consecutive 1 s light pulses (inter–pulse interval of 3 s), with irradiances of 1.6 mW/mm² and a wavelength of 470 nm at two different holding potentials (–40 mV, –80 mV). The experimental data are illustrated in Figure 11A, and the simulation output is in Figure 11B. Once again, they show a close match. This was reaffirmed by comparing these results to the outcomes of the referenced study’s model, which validated the simulation’s accuracy in replicating the experimental data.

6. Discussion

This work presents a simulation framework for a Bio–Neural Dust System (BNDS), providing a valuable tool for advancing research on non–invasive optogenetics. By focusing on computational modeling, this study addresses a critical gap in the experimental validation of neural dust systems, which remains challenging, the resource and time consumption due to its nature, technical issues (e.g., expensive specialized equipment and facility), and ethical constraints. The framework models complex processes involving the piezoelectric M13 bacteriophage, the smooth endoplasmic reticulum (SER), and Aequorin–based bioluminescence. These simulations enable parameter optimization and system behavior prediction, providing a valuable tool for researchers. Furthermore, the software has been made publicly available, allowing others to verify our results and expand on the innovative technology developed in this work.

Our approach builds on well–established methodologies in optogenetics [19] and bionanotechnology [21,22]. Previous studies have demonstrated the feasibility of neural dust systems and the potential of piezoelectric materials for bio–electronic applications [5,6,30]. However, these efforts have primarily focused on hardware development, with limited exploration of integrated system behavior in a virtual environment. For instance, Piech et al. (2020) [6] and Sonmezoglu et al. (2022) [31] developed implantable neural stimulators capable of recording and wirelessly interfacing with neurons via ultrasonic waves, which marked a significant milestone in hardware advancements. The work by Seo et al. (2016) [30] demonstrated wireless recording in the peripheral nervous system using ultrasonic neural dust, highlighting significant progress in implantable neural interface hardware. Similarly, Chorsi et al. (2019) [32] proposed a highly efficient biological sensing device that combines biocompatible materials with piezoelectric components to enable long–term functionality in vivo. These studies primarily emphasize the fabrication, stability, and miniaturization of neural and biological hardware components, which are critical for advancing the field. However, they lack a comprehensive analysis of the integrated system’s behavior when operating under complex biological conditions. There is a notable paucity of comprehensive virtual studies that model the integrated behavior of these systems within complex biological environments. While some computational models exist, they often focus on isolated components rather than the holistic system. For example, Zhang et al.’s (2025) [33] study on piezoelectric transducers provides insights into their individual functionalities but does not encompass the entire neural dust system’s interactions within a biological milieu.

This gap underscores the need for a comprehensive simulation framework, such as the one proposed in this study, which allows researchers to evaluate the integrated behavior of Bio–Neural systems in a virtual environment. By incorporating detailed models for piezoelectric transduction, ion diffusion, and optogenetic responses, our simulator provides a more holistic understanding of these interactions, enabling researchers to optimize system–level performance before transitioning to in vivo experiments.

The simulation framework provides insights into the potential functionality of the BNDS, particularly in non–invasive neural stimulation. It offers a cost–effective and resource–efficient alternative to experimental studies, making it accessible to researchers without specialized facilities. By modifying parameters and integrating additional biological or synthetic components, the simulation framework offers the flexibility to be applied to a wide range of fields, including drug delivery systems and advanced bio–optical sensing technologies. This work opens new possibilities for applications in neuroscience, such as targeted neural modulation and brain–computer interface technologies. The Bio–Neural Dust System (BNDS) offers several advantages compared to other minimally invasive neural interface techniques, such as microelectrode arrays (MEAs), micro–electrocorticography (microECoG), and stentrode arrays. While MEAs provide high spatial resolution and long–term recording capabilities, they are associated with risks of tissue damage and immune responses due to their rigid structures, limiting their biocompatibility and long–term stability [34,35,36]. Similarly, microECoG arrays enable less invasive cortical surface interfacing but are limited by their sensitivity to large–scale neural signals and their lack of depth specificity [37,38,39]. Stentrode arrays, while innovative in utilizing vascular routes for neural interfacing, face challenges in signal resolution and require highly specialized implantation techniques [40].

In contrast, the BNDS harnesses the unique advantages of optogenetics, including high temporal precision and cell–type specificity, combined with neural dust’s minimally invasive nature and biocompatibility. These features make the BNDS particularly promising for applications requiring precise and targeted neural modulation, such as in brain–computer interface technologies and the treatment of neurological disorders. However, challenges such as achieving efficient light delivery and reliable energy transfer still need to be addressed to fully realize the potential of this approach.

While this work focuses on simulation, it aligns closely with experimentally validated principles, such as the piezoelectric effect of M13 bacteriophages [21,22] and the bioluminescent properties of Aequorin [12]. However, experimental validation of the integrated BNDS remains a future goal. Key challenges include fabricating bio–transceivers with the required precision and ensuring biocompatibility in vivo.

Despite its strengths, the framework has limitations. For example, certain assumptions, such as idealized conditions for calcium diffusion and bioluminescence, may need refinement in future iterations. Additionally, extending the framework to model networks of neurons and incorporating more detailed biological dynamics could enhance its utility. Future work will focus on addressing these limitations and validating the framework through experimental studies.

We identify several potential improvements for future development.

• The 3D script of the bio–transceiver is a simplified representation of the mathematical model due to computational limits, making it infeasible to simulate 0.2 pM or more Ca²⁺spheres. Improvements can be made using adaptive particle resolution or Smoothed Particle Hydrodynamics (SPHs), which would enhance realism while managing computational load efficiently.
• The current setup for simulating multiple light pulses necessitates manual code modifications. Future versions should include this functionality within the user interface.
• Expand the neuron model to support networks of neurons and their collective responses to optical stimuli: (a) Introduce spatial coordinates for neurons to model the physical layout of the network. (b) Simulate the effects of light scattering and absorption in tissue, influencing the optical stimuli reaching different neurons. (c) Use light distribution models (e.g., Monte Carlo simulations) to simulate how optical stimuli propagate through the network. (d) Incorporate heterogeneous responses to light based on neuron–specific properties, such as opsin expression levels or membrane dynamics.
• Introduce functionality for exporting simulation data in formats such as CSV for quantitative analysis, JSON for integration with other tools, and PNG or MP4 for visual presentations. These formats enable detailed data analysis, seamless tool interoperability, and clear communication of results in research or presentations.

7. Conclusions

The design and development of this simulation framework necessitated a multidisciplinary approach. Throughout the design and execution stages, we have delineated a high–level architectural plan for the Bio–Neural Dust System and have delivered a foundational version of the simulation framework. Designed with modularity in mind, each component of the framework can be independently altered or enhanced without risking the integrity of the overall system. This paper lays down a robust groundwork that can serve as a springboard for further research and refinement by the scientific community. In closing, our contribution to the simulation framework represents a collaborative stride towards innovation in optogenetics and neuroscience. This work will contribute to the future evolution of non–intrusive brain–computer interface technologies.