Spheroidal Molecular Communication via Diffusion: Signaling Between Homogeneous Cell Aggregates

Recent molecular communication (MC) research has integrated more detailed computational models to capture the dynamics of practical biophysical systems. This paper focuses on developing realistic models for MC transceivers inspired by spheroids – three-dimensional cell aggregates commonly used in organ-on-chip experimental systems. Potential applications that can be used or modeled with spheroids include nutrient transport in organ-on-chip systems, the release of biomarkers or reception of drug molecules by cancerous tumor sites, or transceiver nanomachines participating in information exchange. In this paper, a simple diffusive MC system is considered where a spheroidal transmitter and spheroidal receiver are in an unbounded fluid environment. These spheroidal antennas are modeled as porous media for diffusive signaling molecules, then their boundary conditions and effective diffusion coefficients are characterized. Furthermore, for either a point source or spheroidal transmitter, the Green’s function for concentration (GFC) outside and inside the receiving spheroid is analytically derived and formulated in terms of an infinite series and confirmed with a particle-based simulator (PBS). The provided GFCs enable computation of the transmitted and received signals in the proposed spheroidal communication system. This study shows that the porous structure of the receiving spheroid amplifies diffusion signals but also disperses them, thus there is a trade-off between porosity and information transmission rate. Furthermore, the results reveal that the porous arrangement of the transmitting spheroid not only disperses the received signal but also attenuates it in comparison to a point source transmitter. System performance is also evaluated in terms of the bit error rate (BER). Decreasing the porosity of the receiving spheroid is shown to enhance the system performance. Conversely, reducing the porosity of the transmitting spheroid can adversely affect system performance.


I. INTRODUCTION
Molecular communication (MC) is a bio-inspired mechanism that is envisioned to realize micro-and nano-scale communication systems using molecules as information carriers [1].Despite many efforts by the MC community to model the components of MC systems, more realistic models are required.Existing literature over-simplifies MC components in biological environments and does not sufficiently account for how signaling molecules interact with biological or biosynthetic transmitters and receivers (i.e., cells) or with the environment.Elements and structures of in vitro environments such as Organs-on-Chip (OoCs) could be used to improve model realism, to provide mechanistic insight into the OoC biology, and more generally to contribute to MC research and development.For example, the transmitter could be a group of beta-cells emitting insulin and the receiver could be a group of liver cells that detect the insulin signal and react by increasing its glucose uptake.
Previous MC works have theoretically modeled and simulated the transmitter as a single cell (or machine) releasing molecules with different release mechanisms and the receiver as a single cell detecting the molecules.However, cells do not normally live in isolation but in diverse populations with other cells.This is true in vivo but also in many in vitro systems.In particular, tissues and tumors in multi-cellular organisms and biofilms of microorganisms are common natural instances whereas spheroids, organoids, tumoroids, and islet microtissues are well-known instances in biological experimental setups.Thus, we are inspired by the design of MC transceivers based on a population of (biological or biosynthetic) cells.One realistic transceiver for MC is a spheroid structure, which is a 3D cell aggregation of hundreds or thousands of cells in a spherical shape that is widely used in Organ-on-Chip systems, e.g., in [2].These microphysiological systems have the promising capability to emulate natural organ-to-organ communication dynamics in vitro.For example, a liver-pancreas OoC model with organ-to-organ communication was demonstrated in [2], with the purpose of providing a distinctive platform for studying type 2 diabetes mellitus.In this system, glucose is controlled through communication between liver spheroids and pancreatic islets.The pancreatic islets release insulin under high glucose conditions, which increases glucose uptake by the liver spheroids and reduces their glucose release.The liverpancreas OoC inspires biosynthetic MC systems for healthcare applications.As a first attempt to realize this spheroid-based MC system and study its communication dynamics, we require a suitable biophysical model for an individual spheroid.
Several studies have attempted to model the process of mass transport in a spheroidal environment.The authors of [3] modeled the penetration and diffusion of dyes within multicellular spheroids using diffusion-reaction equations.They assumed constant concentrations of molecules at the inside and outside of the outer spheroid's boundary.Paper [4] focused on mathematically modeling the spatiotemporal dynamics of drugs in spheroids.The authors investigated how drug characteristics impact the propagation of the drug inside a spheroid.As boundary conditions at the spheroid border, they assumed flow continuity and that the concentrations on either side of the outer spherical boundary were equal.The authors of [5] developed an agent-based mathematical model to simulate the passive distribution of microbeads into tumor spheroids while considering the spheroid's growth through externally-supplied oxygen.In their approach, they assumed a constant oxygen concentration at both sides of the spheroid's outer boundary.Paper [6] proposed a model to investigate how factors including the association and dissociation rates, degradation rate, and plasma clearance rate influence the depth of antibody penetration in tumor spheroids.As a boundary condition, they considered the concentration of unbound antibodies inside the spheroid to be equal to that outside multiplied by the porosity parameter or the fraction of tumor volume accessible to the antibody.The authors of [7] developed a mathematical model to examine the diffusion of nanoparticles into tumor spheroids that considers the structural nonuniformity of the spheroid in the radial direction.They applied a boundary condition that determines the concentration inside the spheroid by multiplying the concentration outside with the radiallydependent spheroid porosity.In [8], the researchers presented a mathematical model to predict oxygen diffusion and consumption in tumor spheroids.They employed an analytical solution based on the spherical reaction-diffusion equation with a constant concentration at the tumor boundary that is not dependent on the surrounding concentration.
Furthermore, in our recent study in [9], we mathematically proved that the amplification is indeed present at the boundary and derived the corresponding amplification factor.The provided theorem was investigated through experiments using liver spheroids and glucose as the target molecules.Interestingly, a rapid decrease in the concentration of molecules in the surrounding medium was observed when a spheroid was introduced to a culture medium with a higher glucose concentration.The distinct propagation characteristics prompt us to demonstrate whether spheroid properties facilitate communication.In this direction, we next review relevant existing literature on modeling MC transceiver components.
Transmitters: The most common transmitter model is an ideal point source that releases molecules instantaneously, disregarding physical geometry and realistic release mechanisms [10].Paper [11] pioneered the concept of a pulse-shaped release from a point source transmitter, in which multiple molecules can be released during the pulse.Alternatively, the volume transmitter model proposed in [12] incorporated the transmitter's geometric shape but did not consider a distinct transmitter boundary, leading to an equal concentration inside and outside the transmitter's boundary.Paper [13] studied a box-like transmitter with a surface outlet for controlling molecule release to establish desired concentration gradients.Furthermore, [14] assumed a spherical structure with nanopores for molecule passage.The model in [15] incorporated ion channels on a cell surface transmitter sensitive to either electrical voltage or ligand concentration variations for molecule release.Storage and production dynamics inside the transmitter have also been explored in the literature, ranging from instantaneous [16] to pulse function release [17] and (inspired by neuronal behavior) exponential release based on the number of molecules stored, as proposed in [18].Paper [19] introduced a transmitter model that employs membrane fusion.This mechanism involves the fusion of a vesicle, loaded with molecules and formed within the transmitter, with the transmitter membrane, which enables molecules to be released into the extracellular environment.Moreover, the authors of [20] investigated the controlled release of signaling molecules from compact, non-porous matrix-type drug carriers.The authors also studied the gradual release of molecules from the matrix structure, considering a perfect sink condition at the matrix boundary.However, one can extend their study to porous matrices by incorporating modified diffusion coefficients and drug solubility [21]- [23].
Receivers: The passive receiver, in which information molecules freely diffuse in the receiver's space and the movement of molecules is not affected, is the most common model used in previous works [10].Such a simple model is often used to facilitate analysis of other aspects of an MC system, e.g., the environment boundary [24].The passive model may be relevant for small and hydrophobic molecules (which are repelled from water molecules) that easily pass through a cell membrane.However, many extracellular molecules are too large or too lipophobic to traverse the cell membrane [25] and also cells usually react with signaling molecules (either directly on the surface or in the intracellular environment).
To address the limitations of the passive receiver model, some works have considered a reception mechanism across the receiver (cell) membrane to activate an internal signaling pathway (i.e., a series of chemical reactions controlling a cell function).These works, including [26]- [31], have studied the effects of various reaction mechanisms across the membrane, including cell membrane receptors that can vary in size, number, and spatial distribution.
MC system analysis for microfluidic systems has not investigated relatively large transceivers with different diffusion coefficients.This paper is the initial work in the field to consider spheroidal transceivers under corresponding boundary conditions.However, to establish a complete microfluidic system with spheroids as transceivers, we need to characterize the diffusion and convection of molecules through microfluidic channels and wells, as well as their diffusion and convection within the porous spheroidal structures.In this study, we are primarily concerned with a simplified version of the problem at hand.Hence, in this work, our focus is to model communication between two spheroidal structures, inspired by the cross-talk phenomenon studied in [2].We consider this communication in an unbounded environment without fluid flow.Nevertheless, this assumption is reasonable because the dimensions of the spheroids and the distance between them are significantly smaller than the size of the well in a typical microfluidic system (e.g., 5 mm diameter according to [32]).Moreover, the average flow rate is extremely low in a reference study on which we based our model (i.e., 4.94 µL/min with a complete media turnover time of approximately 2 h according to [2]), supporting our decision to not include convection in this initial work.This research highlights how the distinct features of spheroid structure support the potential application of spheroidal antennas in synthetic biology for transmitting and receiving information.The spheroids' ability to release molecules in a controlled manner enables them to function as transmitters.The porous spheroidal structure also increases the received power of the diffusion signal.Use cases can include improved communication between drug carriers formed from porous matrices for precisely locating diseased cells.
In this paper, we introduce an end-to-end diffusive MC system using spheroids as both transmitter and receiver.The contributions of this paper are summarized as follows: • We propose a novel spheroid-to-spheroid (S2S) communication system in an unbounded environment.• We formulate the molecule release function from the border of a transmitting multicell spheroid as a boundary value problem and derive the corresponding Green's function for concentration (GFC).• As we did for the impulsive point source transmitter in [33], the joint communication channel and spheroidal receiver's response are also modeled as a boundary value problem that accounts for the amplification at the boundary.
• We validate our results with particle-based simulations (PBS), confirming the accuracy and reliability of our model.• We study the geometric effects of the transmitting and receiving spheroids on the molecule concentration inside the receiving spheroid.• We evaluate the system's bit error rate (BER) in the presence of inter-symbol interference (ISI) with different time slot durations and spheroid porosity levels.The rest of the paper is organized as follows.Section II describes the spheroid structure and diffusive MC system using the spheroid.In Section III, the spheroid Green's function is provided.The details regarding the characterization of the diffusive MC channel are presented in Section IV. Results and discussions are outlined in Section V. Finally, Section VI concludes the paper.

A. Diffusive MC System with Spheroidal Transceivers
An end-to-end diffusive MC system is considered, where both the transmitter and receiver are spheroids.We consider two spheroids at fixed locations within an unbounded fluid Fig. 1.System model schematic where the red and greenish-black spheres represent the signaling molecules and cells, respectively.The receiving spheroid is centrally positioned at the coordinate system's origin, and the transmitting spheroid is centered at a distance of rtx along the coordinate system.medium, as shown in Fig. 1.The spherical coordinate system, characterized by (r, θ, φ), is used to illustrate the geometry of the environment, where r, θ, and φ represent the radial, elevation, and azimuthal coordinates, respectively.The origin of the spherical coordinate system is set at the center of the receiving spheroid, while the transmitting spheroid is centered at the arbitrary point rtx = (r tx , θ tx , φ tx ).The transmitter uses signaling molecules of type A. The released molecules move randomly within the environment of the spheroidal transmitter by following Brownian motion before eventually diffusing into the surrounding medium where they continue to freely diffuse.The molecule movements are assumed to be independent of each other.Some of the diffusing molecules will pass into (and potentially out of) the receiving spheroid.
The molecules may bind to cell receptors through a chain of reactions leading to product molecules E that are detected by cells.We do not consider this mechanism in the model and approximate the process based on first-order degradation, i.e., We focus on investigating the concentration profile of type A molecules outside and within the spheroid.
In practice, spheroids are constructed by various methods with the aim to emulate the internal physiological activities of an organ [34].This compact spheroidal structure forms by a spontaneous aggregation of cells where the integrins on their surfaces bind to the extracellular matrix (ECM) of other nearby cells [35].When spheroids reach maturity, they are collected and placed into designated wells.If they grow too large, then oxygen diffusion to the core becomes restricted and this can lead to cell death and the formation of a necrotic center [36].However, in this research, our focus is on small, mature spheroids with a dense arrangement of cells.This means that a signaling molecule moving within the extracellular space (ECS) would always be very close to the cells' boundaries and frequently exposed to bind to receptors.The ECS is the space between cells and is for exchanging nutrients, gases, and signaling molecules between the cells and the surrounding environment.Here, we assume that molecules within the ECS of the spheroid are homogeneously available for transmembrane reactions.We also note that this study does not incorporate the growth and death of cells within the spheroid due to the considerably lower rates of cell growth and death compared to the information transmission rate under consideration [37].

B. Spheroidal Model
As shown in Fig. 1, we consider a spheroid as a threedimensional arrangement of N c cells in a sphere-like shape with radius γ s that mimics the structural and functional characteristics observed in the biological organs or tissue constructs that are used in laboratory environments.The interior space of the spheroid is composed of both the cells and the ECS.The ECS allows for the diffusion of molecules and the formation of concentration gradients.In this paper, we assume that the ECM materials do not contribute a meaningful volume to the occupancy of the spheroid and, consequently, simplify the ECS as empty fluid.Given that V c represents the volume of an individual cell within the spheroid, the volumes of the cell matrix and the extracellular space between cells are equivalent to V c N c and V s −V c N c , respectively, where the overall volume of the spheroid is V s = 4 3 πγ 3 s .In this study, we characterize the architecture of the spheroid as a porous medium, wherein the porosity parameter, ε, serves as the ratio of the extracellular space to the overall volume of the spheroid, i.e., We assume that the spheroids are immersed in a fluid medium that occupies their extracellular space.Signaling molecules, nutrient molecules, and waste products diffuse within the ECS via curved paths, resulting in a shorter net molecule displacement over a given time interval than outside the spheroid.Consequently, the effective diffusion within the whole spheroid volume is reduced compared to the free fluid diffusion outside.We consider the entire spheroid volume to be a homogenized diffusion environment that we will refer to as the "homogenized spheroid" or simply "spheroid" throughout the remainder of this paper.The effective diffusion coefficient within this homogenized spheroid can be determined in terms of the diffusion coefficient (D) of the signaling A molecules in the free fluid as [21] where τ is tortuosity, which refers to the degree of path irregularity or curvature experienced by a molecule while it traverses through the extracellular space of the spheroid, and is modeled as a function of spheroid porosity, i.e., τ =1 ε 0.5 [38].
At the interface between the two diffusive environments with different diffusion coefficients, a continuity condition for flow must be satisfied, which is expressed as and another boundary condition that is generally modeled as [39, Ch. 3] where r ∈ ∂Ω, ∂Ω denotes the spheroid boundary region, Ω is the spheroid region, and c s and c o denote the concentration function inside and outside the spheroid, respectively 1 .The constant κ is a function of the porosity of the medium and should be determined experimentally.In our results in Section V, we demonstrate using a PBS that for two ideal diffusion environments with diffusion coefficients D and D eff , we have κ = D D eff , which was theoretically verified and experimentally studied in [9].Thus, for κ ̸ = 1, a concentration discontinuity (i.e., jump) occurs at the boundary.
Fig. 2 demonstrates the porosity (ε) and boundary concentration ratio (κ) as a function of the number of cells in the spheroid, 15000 < N c < 25000, when the spheroid radius and cell volume are assumed to be R s = 275 µm and V c = 3.14×10 −15 m 3 , respectively [2].As observed in Fig. 2, the boundary concentration ratio increases exponentially with an increase in N c .For N c = 24000, which is the approximate number of cells in HepaRG spheroids [2], we have ε = 0.1349 and correspondingly κ = 4.4919.This value of κ suggests a large concentration discontinuity at the spheroid boundary.

III. GREEN'S FUNCTION FOR THE BOUNDARY VALUE DIFFUSION PROBLEM
To analyze the diffusive spheroidal MC system defined in Section II, we divide the problem into two distinct parts: (i) modeling the transmitted signal given the release of molecules by the spheroidal transmitter with radius γ tx and the effective diffusion coefficient of D tx eff ; and (ii) modeling the received signal at the spheroidal receiver with radius γ rx and the effective diffusion coefficient of D rx eff given the transmitted signal.The Green's function for the boundary value diffusion problem is formulated for each section to characterize the expected diffusion signal.Subsequently, these individual problems are combined into a single problem to analyze the end-to-end diffusive MC system.

A. Diffusion Signal Generated by Spheroidal Transmitter
The spheroidal transmitter is comprised of an aggregation of N tx c cells that are randomly distributed throughout the spheroidal environment to form a homogeneous porous medium.We assume that all cells release molecules simultaneously upon exposure to an external stimulus (e.g., light or sudden change in pH-level).
We begin our analysis by formulating the transmitting spheroid's molecule release rate into an empty environment, i.e., without the receiving spheroid.To simplify the characterization of the diffusion signal from the transmitting spheroid, we assume that the center of the transmitter corresponds to the origin of our coordinate system.First, we obtain the concentration of molecules due to instantaneous release from a point source transmitter positioned at an arbitrary location r0 = (r 0 , θ 0 , φ 0 ) within the transmitting spheroid where one molecule is released at time t 0 .This impulsive point source is represented by the function S(r, t, r0 , t 0 ) = δ(r−r0)δ(θ−θ0)δ(φ−φ0)δ(t−t0) r2 sin θ s −1 m −3 .Given the source S(r, t, r0 , t 0 ), the molecule diffusion inside the transmitting spheroid is described by the partial differential equation (PDE) [40] where c tx (r, t|r 0 , t 0 ) denotes the molecule concentration at point r inside the transmitting spheroid at time t, given the aforementioned impulsive point source.In the spherical coordinate system, the Fourier transform of ( 6) is obtained as where C tx is the Fourier transform of c tx , i is the imaginary unit, and ω is the angular frequency variable in the Fourier transform 2 .The free-space diffusion outside the transmitting spheroid is modeled as In this equation, C otx is the Fourier transform of c otx , which is the molecule concentration at point r outside the transmitting spheroid at time t.The concentration functions c tx (r, t|r 0 , t 0 ) and c otx (r, t|r 0 , t 0 ) that satisfy (7) and (8) subject to the boundary conditions ( 4) and ( 5) are called the Green's functions for concentration (GFCs) of diffusion inside the transmitting spheroid and in the medium, respectively.We have solved the problem and obtained the series-form expression (22) in Appendix A. Due to the linearity of the PDEs, we account for the molecule concentration resulting from all transmitting cells by calculating the integral of c q , q ∈ {tx, otx}, across the transmitter volume Ω TX as where 1 vTX is for normalization.This numerical integral is evaluated with a discrete approximation that improves by increasing the resolution of the discretization.The total summation of molecules inside (N in (t)) and outside (N out (t)) the transmitting spheroid is constant, represented as N tx c × N , where N is the number of molecules released per cell.By normalizing to N tx c × N , we have N in (t) + N out (t) = 1.The release rate is defined as either the negative of the change in the number of molecules inside the transmitter − dNin(t) dt or the change in the number of molecules outside the transmitter . N out is obtained by integrating the concentration of total molecules outside the spheroid.Consequently, the expected release rate g(t) is computed as For tractability in the receiver analysis, we assume that the distance separating the transmitter and receiver is much larger than the sizes of the spheroids themselves.This approximation is consistent with the findings in [12], where the authors compare a spherical transmitter to a point source transmitter and illustrate the improved approximation at greater distances.With this approximation, we can assume that the transmitting spheroid can be represented as a point source with the release function g(t).However, to investigate signal propagation between two nearby spheroids, it becomes essential to consider a volume transmitter with release function g(t) at the spheroid's surface.The investigation of non-uniform, gradient-dependent releases of molecules from individual cells within the spheroid is left for future work.

B. Received Diffusion Signal by Spheroidal Receiver
In this subsection, we formulate the Green's functions for boundary value diffusion problems inside and outside the spheroidal receiver in the absence of the spheroidal transmitter.We assume that an impulsive point source transmitter is positioned outside the spheroid at r0 = (r 0 , θ 0 , φ 0 ), represented by S(r, t, r0 , t 0 ).The free-space diffusion outside the receiving spheroid is described by the partial differential equation (PDE) [40] D∇ 2 c orx (r, t|r 0 , t 0 ) + S(r, t, r0 , t 0 ) = ∂c orx (r, t|r 0 , t 0 ) ∂t . ( In the spherical coordinate system, the Fourier transform of ( 11) is re-written as According to the reaction in (1), the effective diffusion inside the receiving spheroid is modeled as where K is the net A molecule reaction rate due to the binding process characterized by the reaction (1).The boundary conditions at the border of the spheroid are given by ( 4) and (5).We have solved the problem and obtained the series-form expression (22) in Appendix A.

IV. S2S COMMUNICATION
In this section, we first propose the end-to-end diffusive MC channel by determining the molecule concentration received by the spheroidal receiver located at the origin, given the molecules released from the spheroidal transmitter centered at rtx = (r tx , θ tx , φ tx ).Then, we characterize information transmission between the spheroidal transmitter and the spheroidal receiver in an unbounded environment for an on-off keying modulation scheme to send a sequence of binary data.

A. S2S Channel Characterization
We assume that all cells within the spheroidal transmitter release molecules instantaneously at time 0 (i.e., at the beginning of the time slot).This instantaneous release leads to the gradual release of molecules at the boundary of the spheroid, which we approximated with the release function g(t) in (10) from the center of the spheroid.Also, we consider ideal aggregate detection across the entire spheroid, where the observations of individual cells are aggregated and used for spheroid signal processing.Finally, we assume that the mutual impact of the spheroidal transmitter and receiver volumes on the individual GFCs is negligible.
Given c rx (r, t|r tx , t 0 ) in ( 22) as the concentration of molecules at the receiver in response to the impulsive point source located at rtx , the concentration profile given the release function g(t) in the spheroidal receiver, c Total rx (r, t), is obtained as where * denotes the convolution operator.Using the GFC in (14), the probability of observing a specific molecule inside a spheroidal receiver at time t is where Ω RX represents the set comprising all points within the receiving spheroid.The release rate from the transmitting spheroid over ∆t is modeled as a Poisson random variable (RV) with mean (g(t)∆t).
Following [15], we define Y(t) as the number of molecules that are released from the transmitter with release function g(t) within the interval t ∈ [0, T s ], where T s is the time slot duration, and subsequently observed at the receiver at time t within the same interval.Y(t) can be modeled using a Poisson distribution with mean

B. Received Signal in On-Off keying modulation
We adopt the on-off keying modulation scheme where time is divided into time slots of duration T s .In each time slot, the transmission of equiprobable bits 0 and 1 is represented by the average release of 0 and N molecules from each cell at the beginning of the time slot, respectively.At sampling time t s in each time slot, the received signal is defined as the number of molecules of type A inside the spheroidal receiver volume.The sampling time is chosen to maximize the observation probability p obs (t).The receiver utilizes the observed sample to make a decision regarding the transmitted bit.
We define the sequence of transmitted bits for current and W previous time slots, B w = b w ∈ [0, 1], where w ∈ [0 : W ]. When w = 0, it corresponds to the current time slot [0, T s ].Conversely, when w > 0, it indicates the previous time slot [−wT s , −(w − 1)T s ].We immediately have that the number of molecules released at the current time slot (w = 0) and observed by the receiver at time t ∈ [0, T s ] is a Poissondistributed RV with mean b 0 y(t) where y(t) is given by (16).The number of residual molecules released in the previous time slot, w, and observed in the current time slot, w = 0, at time t is denoted as Poisson RV I w (t) with mean The number of molecules observed at time t ∈ [0, T s ], originating from W previously-transmitted symbols, is given by Poisson RV I(t) with mean Therefore, during the current time slot, the receiver's observation at the specific sampling time t s is given by Y R = Y(t s ) + I(t s ).In this equation, Y R is a RV that follows a Poisson distribution with mean b 0 y(t s ) + I(t s ).
The performance of our proposed communication system is characterized by its bit error rate.The bit error rate in (50) is derived in Appendix B.

V. SIMULATION AND NUMERICAL RESULTS
In this section, using the parameters provided in Table I, we first evaluate the boundary concentration ratio, κ, with PBS.We verify the proposed analysis of the transmitting and receiving spheroid Green's functions with this PBS.Additionally, we investigate how the porosity of the transmitting and receiving spheroids affects the molecule concentration inside the receiving spheroid.We compare the received signal at the spheroidal receiver with of a transparent receiver, revealing the signal amplification and dispersion phenomena within the spheroid.Furthermore, we compare the received signal released from the spheroidal transmitter and the signal released from a point source transmitter.Finally, the performance of the S2S communication system is evaluated using the BER.

A. System Parameters and Simulation Methods
The geometric parameters considered in Table 1 are based on the real HepaRG spheroids used in [2] where the spheroid's radius is 275 µm with N c ≈ 24000 cells.The volume of one HepaRG cell is estimated to be 3.14 × 10 −15 m 3 .The molecules reaction rate is set to be k f = 0.01 s −1 .Therefore, the E molecule generation rate inside the spheroid is obtained as ∂c E ∂t = k f c s and correspondingly the term K(ω) in ( 13) is simply k f .
The PBS is implemented in MATLAB (R2021b; The Math-Works, Natick, MA) where the time is divided into time steps of ∆t = 0.5 s, which is sufficiently small for the results to be insensitive to the time step value.The molecules released from the transmitter move independently in the 3-dimensional space.The displacement of a molecule over ∆t is modeled as a Gaussian random variable with zero mean and variance 2D∆t in each Cartesian dimension.For the PBS, the spheroid is simply a diffusion environment with the effective coefficient D eff .When the displacement vector calculated for a molecule outside of the spheroid passes the boundary of the spheroid in a given time slot, the length of the portion of the vector inside the spheroid is updated based on the effective diffusion coefficient D eff , i.e., that part of the vector is scaled by the factor D eff D .Conversely, if the calculated displacement vector of a molecule inside the spheroid in a time slot passes the spheroid boundary, then the length of the portion of the vector outside the spheroid is updated based on the free diffusion coefficient D, i.e., that part of the vector is scaled by the factor D D eff [9].Due to the short time interval (∆t) and the considerable size of the spheroid, the likelihood of a molecule outside the spheroid entering and leaving the spheroid within a single time interval is low.Additionally, considering the substantial distance between the transmitter and receiver, the direct movement of molecules between spheroids within one time interval is improbable.Nevertheless, we still account for these instances of molecules crossing multiple boundaries within one time interval, where we adjust a molecule's displacement vector at every crossing.Considering the reaction (1) inside the receiving spheroid, a molecule may be absorbed by a cell from the extracellular space of the spheroid during a time step ∆t s with approximate probability k f ∆t [41].To represent the homogeneous environment of the spheroidal transmitter, we randomly place N tx c cells within the spheroid according to a uniform distribution.To simulate the release from the transmitting spheroid, we assume that each cell within the spheroid releases N molecules simultaneously in response to a stimulus.To measure a point concentration at a given time and location, we place a transparent sphere centered at that location with a small radius of 10 µm and count the number of molecules inside the sphere at that time.The concentration would then be the counted number of molecules divided by the volume of the sphere.The concentration value is normalized, i.e., is divided by N × N tx c , to obtain the response from the spheroidal transmitter.

B. Single Spheroid Results
Fig. 3 shows the concentration at the inside and outside of the outer spheroid boundary, i.e., c rx (275 − µm, π/2, 0, t) and c orx (275 + µm, π/2, 0, t), respectively, obtained via PBS given the impulsive source at the transmitter and the porosity determined by N rx c = 24000 cells within the spheroid.We have also plotted the concentration at the outside of the boundary scaled by the factor κ to verify the boundary condition (5).The minor mismatch at the peaks is mainly due to the PBS procedure to approximate the concentration at a point.To compute the concentration at the boundary point (275 µm, π/2, 0), we have assumed a transparent sphere of radius 10 µm centered at the boundary point.By partitioning the sphere into left and right hemispheres, we have counted the molecules within each hemisphere.Subsequently, we have normalized the counts by dividing them by the hemisphere volume to approximate the concentration of molecules at the inside and outside of the boundary.The curvature of the boundary causes slight variations in hemisphere volumes, which leads to slightly higher and lower concentration approximations outside and inside the boundary, respectively.
Fig. 4 depicts the concentration at different points inside the spheroidal receiver given the instantaneous point transmitter at r0 = (500 µm, π/2, 0) obtained from the analysis provided in Appendix A in (22) and also the PBS when N rx c = 24000.As observed, the PBS fully confirms the analytical results.The figure also indicates that the concentration signal weakens significantly at the points closer to the spheroid center and farther from the point source.In this case, we observe that the peak concentration at the center is about 8 times smaller than that at the boundary.
According to the findings in Fig. 3 and Fig. 4, as molecules penetrate the porous structure of the spheroid, they become trapped in the pores near the border, which restricts (but does not prevent) their movement inward or outward.This intrinsic spheroid property proves advantageous for MC.The receiver's spheroidal structure enables it to absorb molecules from the environment and keep them within its porous structure for an extended period, which leads to an amplified diffusive signal.

C. End-to-End S2S Results
Fig. 5 shows the molecule release rate from the transmitting spheroid with N tx c = 24000 cells as obtained from (10).As time increases, the release rate decreases since the dynamics of molecules diffusing out of or back into the transmitter are more extreme at the beginning.The analytical release rate was also consistent with the PBS.
In Fig. 6, we illustrate the concentration profiles of molecules released from the spheroidal transmitter with release function g(t) and observed in the center of the receiving spheroid with and without degradation as obtained from (14).As observed, our analytical results for the communication between two spheroids are verified with PBS.By increasing the degradation rate inside the receiving spheroid, we have a lower concentration.For clarity of exposition, we only consider analytical results for the remaining figures in this paper.
We have provided Figs.7 and 8 to explore the impact of spheroid porosity on signal propagation in the S2S system.In Fig. 7, we keep the porosity level constant in the transmitting spheroid and vary the receiving spheroid porosity via (2) by adjusting the number of receiving cells in order to examine how it affects the total concentration of molecules inside the receiver obtained from (14).The results show that reducing the porosity of the receiving spheroid inhibits the rapid molecule diffusion within the porous medium, which leads to a delayed

PBS Analysis
Fig. 6.The molecule concentration released from a transmitting spheroid with porosity εtx = 0.1349 at a distance of 1000 µm, evaluated at the center of the receiving spheroid with porosity εrx = 0.2791, is compared for two different degradation rates, k f , over time.The data is obtained through analysis and PBS simulations for N tx c = 24000.
increase in the signal.This is because, upon penetration into spheroids with lower porosities, the molecules get stuck inside the more tortuous path between cells and experience a lower effective diffusion coefficient.Their limited ability to escape the porous environment leads to signal amplification.As observed, the spheroidal receiver model displays signal amplification compared to the transparent receiver but at the cost of delay, which can impact the transmission rate.In Fig. 8, we maintain a constant porosity level in the receiving spheroid to examine the influence of the porosity of the transmitting spheroid on the total concentration of molecules inside the receiver as obtained from (14).To isolate the impact of the transmitting spheroid's porosity, we normalize the GFCs by dividing it by the total number of molecules released.Hence, we assume constant "power" (i.e., number of molecules released) per transmitter regardless of porosity value.As we increase the porosity level of the transmitter, molecules find it easier to move within the spheroid, causing them to escape more quickly and easily from the porous structure of the Fig. 9. BER of S2S system as a function of Ts for various receiving porosities over a distance of 1000 µm, with εtx = 0.1349.The degradation rates are set as k f = 0 for the transmitting spheroid and k f = 0.01 for the receiving spheroid.Each cell within the spheroid releases only one molecule to transmit a bit 1.We also compare the results with the transparent receiver.Fig. 10.BER for S2S system as a function of the number of cells in the transmitting spheroid over a distance of 1000 µm, with εrx = 0.2791.The degradation rates are set as k f = 0 for the transmitting spheroid and k f = 0.01 for the receiving spheroid.The number of released molecules per transmitter is fixed at 24000 molecules.We also compare the results with the point source transmitter.spheroid.This leads to a reduction in delay and an increase in the peak of the received signal.Additionally, we compare our spheroidal transmitter with an instantaneous point source transmitter, which indicates the lowest delay and the highest peak signal since the molecules do not get stuck in the porous environment of the transmitting spheroid.
We also evaluate the communications performance of our model by plotting the probability of error in terms of time slot duration for different receiving spheroid porosity in Fig. 9. Given the considerable cell population within the transmitting spheroid, we assume that each cell within the spheroid releases only one molecule to transmit a bit 1 to enable a meaningful performance analysis.Our findings reveal that if we have lower spheroid porosity, as illustrated in Fig. 7, then we have an amplified signal in the receiving spheroid, which can decrease the BER even though the signal is delayed.Furthermore, we conducted a BER comparison between the S2S system and a system where the transmitter is a spheroid and the receiver is transparent.The result demonstrates that incorporating the porous structure for the receiver leads to a lower BER compared to using a transparent receiver.The signal amplification caused by the spheroidal structure of the receiver demonstrates effectiveness in improving the accuracy of the detection process within the receiver.Finally, we measure the impact of transmitting spheroid porosity on the BER in Fig. 10 by adjusting the number of transmitting cells.By increasing the number of cells or decreasing the porosity, with a fixed number of released molecules per transmitter, we observe a lower molecule concentration in the receiving spheroid with higher delay, as shown in Fig. 8.This weakened and delayed signal is associated with a subsequent increase in the BER.The figure also compares the BER between the S2S system and a point source transmitter releasing 24000 molecules per bit 1.This analysis underscores that as the porosity of the transmitting spheroid decreases in the S2S system, the system's perfor-mance approaches that of a system comprised of a point source transmitter and a spheroidal receiver.

VI. CONCLUSION
We proposed an end-to-end diffusive molecular communication (MC) system using a spheroidal transmitter and receiver.We successfully derived the analytical solution for the Green's function for concentration (GFC) for both transmitter and receiver spheroids.As confirmed by PBS, we investigated how different porosity values in both the spheroidal transmitter and receiver impact the GFC and consequently the received signal.Our findings indicate that increasing the porosity of the receiving spheroid leads to a lower peak GFC with a reduced delay.In addition, we observed that increasing the porosity of the transmitting spheroid while maintaining a constant number of released molecules from the entire spheroid results in a higher concentration peak inside the receiver and a shorter delay.The evaluation of our proposed system's performance is conducted using the bit error rate (BER) in terms of time slot duration and transmitter porosity under the effect of intersymbol interference (ISI).The results show that lower porosity in the receiving spheroid, although causing signal delay, leads to reduced BER since it amplifies the concentration, which can enhance the detection accuracy in the receiver.However, by decreasing the porosity of the transmitting spheroid while keeping the transmission power constant, the receiving spheroid experiences a decrease in molecule concentration, resulting in a higher BER.
This paper represents an initial step toward an MC system that employs a transmitter and receiver model based on aggregations of cells.The results of this research will create a solid foundation for future advancements in the field of MC.In future work, one might improve model accuracy or examine signal diffusion in a nonhomogeneous porous media model of the spheroid, for example by including (i) the necrotic core that can result from cell death in the center of a spheroid due to a lack of oxygen; or (ii) a spheroid of two (or more) different cell types, e.g., HepaRG cells and stellate cells.Another challenge is to perform experimental validation, which would require measurements of internal spheroid concentration.Future studies can consider extending the model to more than two spheroids to illustrate more realistic scenarios.Considerations can also include modeling spheroidal transceivers within bounded microfluidic wells, accounting for fluid flow dynamics that could promote deeper molecule penetration into the spheroid.

APPENDIX A TRANSMITTER AND RECEIVER ANALYSIS
In the following, we obtain analytical expressions of Green's function for both the spheroidal transmitter and receiver in the frequency domain.We consider homogeneous initial conditions for the system.The Green's function for the boundary value diffusion problem for inside and outside of the spheroidal transmitter is given by ( 7) and ( 8), given the boundary conditions ( 4) and ( 5) by replacing c s and c o with c tx and c otx , respectively.Also, the Green's function for the boundary value diffusion problem for outside and inside of the spheroidal receiver is given by ( 11) and ( 13), given the boundary conditions ( 4) and ( 5) by replacing c s and c o with c rx and c orx , respectively.

A. General Analytical Solution
The general solution of both problems is considered as [42] C q (r, θ, φ, ω|r 0 ) = ∞ n=0 n m=0 H mn R q n (r, ω) cos(m(φ − φ tx )) × P m n (cos θ), q ∈ {tx, otx, orx, rx}, where R q n (r, ω) is the unknown radial Green function of r for the region q ∈ {tx, otx, orx, rx}, F(θ, φ) = cos(mφ)P m n (cos θ) is a spherical harmonic with degree n and order m that satisfies the partial differential equation (PDE) in and H mn denotes the unknown coefficient that needs to be determined for the particular configuration.Also, the representation of δ(φ − φ 0 ) δ(θ−θ0) sin θ based on the aforementioned spherical harmonic is given by where L 0 = 1 2π and L m = 1 π , m ≥ 1.

B. Transmitter Region
Replacing C q and δ(φ − φ 0 ) δ(θ−θ0) sin θ in ( 7) and ( 8) by the corresponding series-form representations given by ( 22) and (24), respectively, leads to (19) and (20) at the top of the next page.Matching the two sides of (19) yields and where where k otx = −iω D .By applying (22) in the Fourier forms of boundary conditions (4) and ( 5), we obtain where is the concentration ratio at the boundary of the transmitting spheroid.In these equations, superscripts + and -represent a minuscule distance above or below the transmitter radius γ tx , respectively.
To establish the discontinuity boundary condition of the Green's function at the point source location r = r 0 , we integrate (26) across a small interval [r − 0 = r 0 −ϵ, r + 0 = r 0 +ϵ].Taking the limit ϵ → 0 and employing the sifting property of the Dirac delta function at the right-hand side leads to the boundary condition [42] r 2 ∂R tx n (r, ω) ∂r The solutions of the homogeneous form of the PDE of (26) and also (27), intended for the transmitter, are then given by where j n (•) and y n (•) are the nth order of the first and second kinds of spherical Bessel function, respectively, and h n (•) represents the nth order Hankel function of the first type.By applying the solutions ( 31) and ( 32) to the boundary conditions ( 28), (29), and (30), and also the continuity condition of concentration at r0 , i.e., R tx n (r, ω) the system of linear equations ( 21) is obtained at the top of this page and from which the coefficients A n , B n , G n , and D n are calculated.We note that j ′ (•) and y ′ (•) in ( 21) are the derivatives of the functions j(•) and y(•) with respect to r, respectively.By having the coefficients A n , B n , C n , and D n , R q n (r, ω), q ∈ {tx, otx} in ( 31) is known and the GFCs of the inside and outside of the transmitting spheroid are computed by taking the inverse Fourier transform of ( 22) for the inside and outside, respectively.

C. Receiver Region
Similar to the transmitter case, replacing C q and δ(φ − φ 0 ) δ(θ−θ0) sin θ in ( 12) and ( 13) by the corresponding series-form representations given by ( 22) and ( 24), respectively, leads to (34) and ( 35) at the top of the next page.Matching the two sides of (34) yields and where k rx = −iω D .Similarly, ( 35) is reduced to where . We emphasize that k orx is a function of the degradation rate via K.
Similar to the transmitter scenario, when we apply (22) in the Fourier forms of boundary conditions (4) and ( 5), above or below the receiver radius γ rx , we obtain R rx n (r, ω) where κ rx = D D rx eff is the concentration ratio at the boundary of the receiving spheroid.
Similar to what we did to obtain (30), we integrate (38) across a small interval below and above the point source y n (k orx r 0 ) j n (k orx r 0 ) −h n (k orx r 0 ) 0 r 2 0 y ′ n (k orx r 0 ) location r = r 0 , which leads to the corresponding boundary condition [42] r 2 ∂R orx n (r, ω) ∂r r=r + 0 − r 2 ∂R orx n (r, ω) ∂r The solutions of the homogeneous form of the PDE of (38) and also (39), applicable to the receiver region, are then given by R rx n (r, ω) = G n j n (kr), r < γ rx , R orx n (r, ω) = B n j n (kr) + A n y n (kr), γ rx < r < r 0 D n h n (kr), r > r 0 . ( By applying the solutions (43) and (44) to the boundary conditions (40), (41), and (42), and also the continuity condition of concentration at r0 , i.e., R orx n (r, ω) the system of linear equations ( 36) is obtained at the top of this page and by which the coefficients A n , B n , G n , and D n are calculated.By having the coefficients A n , B n , G n , and D n , R q n (r, ω), q ∈ {rx, orx} in ( 43) is known and the GFCs of the inside and outside of the receiving spheroid are computed by taking the inverse Fourier transform of (22) for the inside and outside, respectively.

APPENDIX B ON-OFF KEYING S2S PERFORMANCE ANALYSIS
The RV Y R , which follows a Poisson distribution, can be described in the following manner: where notation Pr(•) and E(•) are used to represent the probability function and the expectation operator, respectively.For error probability analysis, we adopt the genie-aided decision feedback (DF) detector [43].This detector benefits from the assistance of a genie, which informs it about the previously-transmitted bits, i. (48) Here, B0 represents the estimated transmitted bit in the current time slot.In practice, the previously-transmitted bits, B w = b w for j ∈ {1, ..., W }, are unknown.As a result, previous decisions, Bw = bw for w ∈ {1, ..., W }, need to be utilized in (46) instead.
We can derive the threshold decision rule for receiver output in the current time slot, y, from (48).According to this rule, if y ≤ ξ, then B0 is set to 0, and if y > ξ, then B0 becomes 1.We calculate ξ as follows [44] To compute this threshold, having perfect knowledge of the channel, i.e., p obs (t s ), is necessary.The error probability of this detector can be determined by the following expression: (52)

Fig. 5 .
Fig. 5. Molecule release rate from the transmitting spheroid for N tx c = 24000.

Fig. 7 .Fig. 8 .
Fig.7.Impact of receiving spheroid porosity on average received concentration over the entire spheroid from a transmitting distance of 1000 µm, where εtx = 0.1349.The degradation rates are k f = 0 for the transmitting and receiving spheroids.We also compare the results with the transparent receiver.