Heuristic Barycenter Modeling of Fully Absorbing Receivers in Diffusive Molecular Communication Channels

In a recent paper it has been shown that to model a diffusive molecular communication (MC) channel with multiple fully absorbing (FA) receivers, these can be interpreted as sources of negative particles from the other receivers’ perspective. The barycenter point is introduced as the best position where to place the negative sources. The barycenter is obtained from the spatial mean of the molecules impinging on the surface of each FA receiver. This paper derives an expression that captures the position of the barycenter in a diffusive MC channel with multiple FA receivers. In this work, a heuristic model inspired by Newton’s law of gravitation is found to describe the barycenter, and the result is compared with particle-based simulation (PBS) data. Since the barycenter depends on the distance between the transmitter and receiver and the observation time, the condition that the barycenter can be assumed to be at the center of the receiver is discussed. This assumption simplifies further modeling of any diffusive MC system containing multiple FA receivers. The resulting position of the barycenter is used in channel models to calculate the cumulative number of absorbed molecules and it has been verified with PBS.

models to calculate the cumulative number of absorbed molecules and it has been verified with PBS data in a variety of scenarios.

Index Terms
Diffusive molecular communication, barycenter, fully absorbing receiver, multiple receivers.

I. INTRODUCTION Molecular Communication (MC) is an interdisciplinary communication paradigm that relies
on molecules propagation to exchange information.This unique discipline opens the door to establishing communication on the scale of nanometers to micrometers that can be used between nanorobots or investigating and controlling the natural communications occurring around us.MC studies will lead to reliable cooperation between nano-devices to increase the complexity of their tasks.MC can be divided into two different classes known as natural and artificial.Natural MC has evolved over millions of years to perform various functions in biological systems and there is an excellent potential to investigate it from the communication and information exchange point of view [1].At the same time, artificial MC is a human-made field that seeks to develop communication systems based on the principles of natural MC.Natural MC occurs in biological systems, where molecules such as hormones, neurotransmitters, and pheromones are used to transmit information between cells, organs, and individuals [2].For example, in the nervous system, neurotransmitters such as dopamine and serotonin are released by neurons to transmit signals to other neurons or to muscle cells [3].In the immune system, cytokines and chemokines are released by cells to signal the presence of pathogens or tissue damage [4].Artificial MC systems, on the other hand, involve the design and implementation of artificial molecules, often using nanotechnology, to transmit information between artificial entities such as sensors, robots, and implants [5].In these systems, the artificial molecules are designed to mimic the behavior of natural molecules and interact with artificial receptors to transmit information.One of the advantages of MC is its potential for use in environments where electromagnetic communication is not possible or desirable.For example, MC can be used in applications such as targeted drug delivery, nanomedicine, and implantable devices, where electromagnetic radiation can be harmful or interfere with the operation of the device [6], [7].Additionally, molecular communication can be used in underwater environments where electromagnetic waves have limited range and are subject to interference [8].MC has been studied under different conditions, e.g. with active or passive receivers, instantaneous or temporal release of molecules, different boundary conditions of the physical channel and so on [9].

A. Previous Works
Almost all the works published so far have analyzed Single Input Single Output (SISO) model.
However MC systems are intrinsically Multiple Input Multiple Output (MIMO).In fact an MC system with multiple receivers is closer to reality since an individual receiver alone can not afford the complexity of the task.In natural MC that are occurring around us, multiple receivers always cooperate as a unique system.Therefore, it would be a considerable achievement to study multiple receivers and consider their interaction.
Modelling the MC with multiple fully absorbing (FA) receivers is a complicated task due to the interaction between the receivers.Some attempts have been made to model the effect of multiple FA receivers in the channel.Authors in [10] tried to model the channel's impulse response using a function similar to the SISO system response by applying curve fitting algorithms.Specifically, the paper's channel modeling section used control coefficients over the SISO model to describe MIMO system.The Control coefficients were selected to comprehend system characteristics and being time independent.Then, nonlinear regression models were used to fit simulation data.Similarly [11] relied on the simulation results under different scenarios, and proposed the empirical formulas of the cumulative absorbing probability and the absorption probability density on the intended receiver with respect to the angle between two receivers, distances from transmitter to the receivers, and the spherical receiver size.However, there is no guarantee that such approach can ensure generality of the resulting model.
Authors in [12] endeavored to derive an analytical model for a Single Input Multiple Output (SIMO) system operating under conditions of negligibly small mutual interaction between the transmitter and receivers.Such a configuration is tantamount to a scenario wherein the distances separating the receivers from the transceiver, as well as between the receivers themselves, are of a sufficiently great magnitude.However, the key missing part of the generic model for the presence of multiple receivers is that it captures the receivers' observation for any arbitrary positions.They did not consider the cases where receivers are blocking the line of sight of one another or when a receiver is close to the source.Because in that case, the system of the equation describing the model needs to be solved by numerical integration due to strong correlations appearing among the receivers.Moreover, the main limitation of that work and all the related papers available in the literature [13], [14] is that they miss paying attention to the interaction among the receivers when they are close to the source.To model the diffusive MC system with multiple FA receivers they can be substituted with the negative point sources and their positions are called the barycenters.Previous works did not take into account where to put the negative source.They just put it in the center of the FA receivers (although this false assumption had been discussed and showed thoroughly in the empirical study of [15] and will be shown with proofs and mathematical discussions in this paper).
The MC channel model with multiple FA receivers has been recently proposed in [15].The authors proposed an analytical model to describe the impulse response of the diffusive channel between a point transmitter and a given number of FA receivers in an MC system.The presence of neighboring FA nanomachines in the environment was taken into account by describing them as point sources of negative molecules.A fundamental problem was the question: "Where should the negative point source be placed?"The authors gave an answer to this dilemma by defining the barycenter point, which is the spatial mean of the molecules that have hit the surface of the FA receiver.The authors in [15] developed an empirical expression to describe the position of the barycenter by applying curve fitting.They found that if there is a source and a receiver in the environment, then the position of the barycenter lies between the center and the surfacepoint of the spherical receiver.The surface-point is defined as the closest point on the surface of the receiver to the transmitter.They proposed a parameter γ as a function of time and distance between the source and the receiver.The proposed parameter varies between one and zero.When γ is equal to one, it means that the spatial average of the molecules from the source are all concentrated at the surface-point of the receiver, while γ equal to zero means that the molecules are distributed around the receiver, and their spatial average coincides with the center of the receiver.To obtain an analytical barycenter, one must know the distribution of the particles on the surface of the receiver.

B. Contributions
In this paper, we derive a heuristic analytical expression that locates the barycenter point of the spherical FA receiver.Having these results allow us to obtain a model that can capture the expected number of absorbed molecules by multiple FA spherical receivers at any arbitrary position.Moreover, the model allows us to have an understanding of under which circumstances we can apply simplifying assumptions and skip the computation of the barycenter and assume it is located at the center of the receivers.This brings the analogy with the common far-field assumption existing in conventional electromagnetic-based communication.First, we describe the system model according to [15], then the derivation of the barycenter is shown.Ultimately, we compare the number of molecules absorbed by the receivers based on the resulting model with the analytical barycenter, the empirical barycenter calculated from the particle-based simulation (PBS), and the cumulative number of absorbed molecules obtained directly from it.We considered a single transmitter and two FA receivers MC scenario to ensure that our contribution can be used to describe the presence of a second receiver around the intended receiver.We also compare the analytical γ proposed in this paper with the empirical γ computed by using the PBS data.Simulation results for the case of two receivers with different radii and five receivers in close proximity are also shown.Hence, despite other available papers in the literature which considered specific scenarios valid for certain conditions, in this paper, we discuss the most complicated scenarios to prove the generality of our model.We believe that using the tools and the methodology introduced in this paper is the missing part to step forward toward modeling the diffusion-based MC with multiple FA receivers.

C. Outline
The rest of the paper is organized as follows.Sec.II describes the system model and introduces the negative point source.Then, it starts from SISO modeling and extends it to the scenario with two FA receivers and finally multiple FA receivers.In Sec.III we propose the analytical barycenter model, verify it, and investigate its behavior.Sec.IV illustrates the simulation results and validates the model.Finally, Sec.V provides the conclusion of this paper.

II. SYSTEM MODEL
In this section we discuss the system model for three cases.To begin with, we review the SISO model.Then, we develop a Single Input Two Output (SITO) model and introduce the concept of negative source.Lastly, we extend the SITO case to a SIMO scenario.The transmitter is pointwise and emits N T messenger molecules of the same type into the environment instantaneously.
The molecules emitted by the transmitter diffuse with constant diffusion coefficient D µm 2 /s through the medium between transmitter and receivers in an unbounded 3D environment.The receivers have FA properties with a spherical geometry and are able to count the number of absorbed molecules.Once the molecules hit the surface of the receiver, they stop moving.The FA property leads to a coupling effect between the receivers.So to study the number of molecules absorbed by the receivers, we have to take into account the interaction between the receivers.

A. Single Input Single Output (SISO)
Diffusive molecules propagation is governed by Fick's second law that links the time derivative of the flux to the Laplacian of the molecules' concentration p (r, t) at distance r and time t as [16] ∂p (r, t) ∂t = D∇ 2 p (r, t) .
The initial and boundary conditions of (1) vary depending on the MC system characterization.
Yilmaz et al. [17] specified the boundary and initial conditions as an impulsive release of molecules, unbounded environment, and an FA spherical receiver R with radius R.They obtained the expression for the hitting rate of the molecules on the surface of the receiver, namely f d (C,T ) , t , which depends on the distance d (C,T ) between the transmitter T and the center of the receiver R, at time t.The channel impulse response of a diffusive MC channel with a single spherical FA receiver of radius R centered at distance d (C,T ) from the transmitter reads and the absorption rate, i.e., the number of molecules absorbed by the receiver per unit time is when the transmitter T emits N T molecules impulsively.The number of absorbed molecules is obtained by integrating (3) up to time t where is the complementary error function.

B. Single Input Two Output (SITO)
When two receivers are present in the same channel, their absorption rate no longer follows (3) due to their full absorption characteristic.The presence of a second FA receiver has the effect of removing molecules from the environment, thus reducing the absorption rate of the first receiver.
The coupling effect of FA receivers on each other has been taken into consideration by introducing the concept of the negative point source of molecules.To model the negative source effect of receivers, we can consider the existence of a negative point source and replace it with the FA receivers, except for the desired receiver, which must be investigated in terms of the number of absorbed molecules.With reference to the receiver R 1 , its hitting rate is influenced both by the number of molecules released by the transmitter and by the reduction of molecules in the environment due to the presence of the receiver R 2 .From the R 1 perspective, R 2 can be interpreted as a point source of "negative" molecules, characterized by the fact that the number of released molecules coincides with the absorbed ones up to a given time.As shown in [15], the best position to place the fictitious negative point source is given by the absorption barycenter point.The barycenter of each receiver is defined as the spatial average of the molecules that adhere to the surface of the receiver due to the FA property.The mutual interaction between the two FA receivers can be modeled by applying the superposition principle.
Mathematically, the hitting rate (2) must be combined with an expression that solves (1) and satisfies the additional boundary condition at R 2 , which absorbs molecules with an (unknown) absorption rate n 2 (t).The effect of this absorption is accounted for as a negative source.The effect of negative source signal, which is concentrated in the absorption point, perturbs n 1 (t) according to the channel impulse response (2).Obviously (2) is the response to an impulsive release.Since n 2 (t) varies with time, we take the convolution between them.Because of the symmetry, we can apply the same reasoning by swapping the roles of the absorption rates n 1 (t) and n 2 (t).We can evaluate the absorption rates of the two receivers using where denotes the convolution, is the distance between the center of R 2 and barycenter of R 1 .
To determine the expected cumulative number of absorbed molecules on receivers, the integral of ( 6) is required.Taking the Laplace transform of the integration of (6), one obtains where L {f } = f and L {N } = N .We can write (7) as a matrix multiplication Thus the solution in the S domain obtained by a matrix inversion followed by multiplication Applying the inverse Laplace transform on (9) results in that expresses the expected cumulative number of absorbed molecules by R 1 , where R 1 and R 2 are the radius of receivers R 1 and R 2 respectively [15, eq. ( 21)].

C. Single Input Multiple Output (SIMO)
Following the same reasoning as for the SITO case, the system of equations corresponding to multiple receivers in the S domain can be written as Applying the Laplace transform on (11) allows us to write it in terms of matrix multiplication Unlike the SITO case, the time domain closed-form solution of ( 12) has not been derived yet.
However, since ( 2) is causal, we can solve the system of equations numerically.

III. THE BARYCENTER ANALYTICAL MODEL
The main objective of the present paper is to obtain an analytical expression that locates the barycenter of all the receivers an MC system with multiple FA receivers.Because knowing the position of the barycenter allows us to replace other FA receivers with negative point sources and consequently solve (11).In the following, we model the barycenter point when there are two FA receivers in the channel by taking advantage of the results from [18] and the superposition principle.Then we extend the model to the case of an arbitrary number of FA receivers in the channel.

A. Barycenter in SITO
The barycenter point is the average of the position of molecules that hit the surface of the receiver up to time t.It is located inside the receiver's volume and depends on time and the position of the transmitter with respect to both the intended and the other receiver.We define it as the weighted sum of two vectors in 3D space such that the first describes the effect of the transmitter and the second the effect of the other receiver on the intended one.Hence the position of the barycenter point inside the volume of receiver R 1 can be written as where the boldness of a symbol indicates that it represents a vector in 3D space, and corresponds to the barycenter point as the effect of the transmitter on R 1 and Let us investigate B (R 1 ,T ) assuming that there is no other receiver around.For a spherical receiver, it is clear that B (R 1 ,T ) varies on the radius of the sphere towards the transmitter.
where the parameter γ specifies that the point B (R 1 ,T ) is located between the surface-point and the center of the spherical receiver R 1 .Derivation of the γ is explained in the following.
To achieve γ we need an expression that describes the distribution of the particles over the surface of the receiver.Authors in [18] modeled the diffusive MC channel with a single spherical receiver in a spherical coordinate assuming that the receiver is located at the center of the coordinate.However, their goal was not to describe the distribution of the particles over the DRAFT April 27, 2023 surface of the receiver.We used their resultant derivation and combined it with our interpretation to have an understanding of the distribution of the particles over the receiver and consequently obtain the γ.They defined the concentration of the molecules as p(r, θ, φ, t) thus, the diffusion equation and its corresponding boundary and initial conditions are The instantaneous release of the molecules from the point source into the environment at time t → 0 is defined as The unboundedness of the environment is represented as The reaction at the surface of the receiver, r = R, can be described by where w is the reaction rate at the surface of the receiver.To characterize a fully absorbing receiver w should tend to infinite.The solution of the diffusion equation in terms of the cumulative number of absorbed molecules was obtained in spherical coordinate as [18, eq. ( 33)].
The solution was written in terms of an integration from 0 to π.This observation has inspired us to create the surface of a sphere by putting an infinite number of rings as shown in Fig. 1 to construct the surface of a spherical receiver.
By looking closely at the integration [18, eq. ( 33)], we can decompose the integral based on the angle θ.Consequently, the fraction of absorbed molecules on rings that create the surface is where But ( 19) depends on w and in order to use it in the case of FA receivers we need to find out Y as w → ∞ (See Appendix A).In the presence of only a transmitter, the spatial average of absorbed particles on the surface of the receiver is expected to be on its radius that is oriented towards the transmitter.This property is due to the symmetrical propagation of particles in all directions around the line of sight.To compute γ, we consider the absorbed particles of each ring on the surface of the receiver as the weight of the points on the diameter that includes the aforementioned radius.Those points are also the center of each ring.We can take the weighted sum of the coordinate of the circles' centers on the diameter of the sphere aligned with the transmitter as The numerator of (24) is the weighted sum of the points on the diameter of the receiver, which includes the radius of the receiver in the direction of the transmitter.The idea of creating the DRAFT April 27, 2023 integral on the numerator is that we want to sum the 1D coordinate of the points i.e., ranging in [−R, R], on the specified diameter such that each point has a weight that comes from a ring on the surface of the sphere.Note that the aforementioned points are the centers of the rings.To this aim, we take the integral on θ in the range of [0, π].The relation between the rings and their centers on the specified diameter can be expressed by R cos (θ), which indicates the center of the rings as shown in Fig. 1.Hence we take the integral of all rings as weights of their centers on the diameter.The denominator acts as a normalizing factor because we want that the value of γ stays between one and zero.This bound is equivalent to the introduced concept before that the barycenter point in a SISO system remains on the radius between the center (γ = 0) and surface-point (γ = 1) of the receiver.
Fig. 2 compares the empirical γ obtained from the PBS with the analytical γ from (24).The PBS data was obtained through a Monte Carlo simulation for 100 trials with simulation step time 10 −7 s and diffusion coefficient 79.4 µm 2 /s.The radius of the receiver is assumed to be 1 µm and we change the distance between the transmitter and the center of the receiver from 1 to 12 µm.The duration of PBS is 2 s.As we expected by increasing the distance between the transmitter and the receiver the value of γ reduces because the signaling molecules are spread in the space and hit the FA receiver more uniformly compared to the scenario with a close distance between the transmitter and receiver.We would like to note that the slight difference between the empirical γ and the analytical one is due to the time discretization of the PBS [19].Although we chose simulation step time equal to 10 −7 s even in this case, not all the particles are trapped over the surface of the FA receiver.Some trap inside the receiver during the PBS.This phenomenon appears due to the discretization of the time domain in PBS while in continuous diffusion we assume particles are absorbed as they hit the surface of the receiver.We would like to underline that the slight observed mismatch is a consequence of temporal discretization in PBS.Finding a solution to deal with the mismatch between Brownian motion and continuous diffusion with absorbing boundary conditions is out of the scope of this paper.After modeling the effect of transmitter on the receiver, we study the effect of R 2 , on R 1 , which was shown by B (R 1 ,R 2 ) .We follow a similar approach as the effect of transmitter since we assume that the other receiver performs as a negative point source from the R 1 's perspective.
However, in this case instead of having an attraction effect, which was locating the barycenter between the center of R 1 and the surface-point, due to the negativity of the source we assume a repulsion effect.Thus the barycenter point as a result of the negative source is somewhere between the center and the farthest point on the surface of the sphere from the center of R 2 .A simpler way to formulate this definition is just changing the sign of γ where is the point on the surface of R 1 towards the center of R 2 .
We need to combine the contribution of the two sources but it must be taken into account that they do not have the same power to define the position of the barycenter.The effect of the two sources (i.e.positive and negative) should not be the same because the number of molecules released by the transmitter is different from that absorbed by the other receiver, which is equivalent to the number of negative molecules.Hence, we consider the coefficients ζ (1,1) and to equalize the contribution from each source.We propose the value of 1 to describe the transmitter's effect on the barycenter.On the other hand, to model the effect of the other FA receiver, which is modeled as the fictitious negative source, we have been inspired by Newton's law of universal gravitation [20].Let us consider each receiver as a planet in the universe.
Hence the planets can have gravitational forces on each other that is inversely proportional to their squared distance.Moreover another important factor is their mass.But in our case a reasonable parameter that can represent the effective mass of the planets (FA receivers) is their radius.Finally we propose the coefficients to take into consideration the interaction between the receivers as ).By normalizing the coefficients we have The barycenters of receivers in a SITO system becomes By looking at ( 14) and ( 25) we understand that the barycenter is highly dependent on the behavior of γ contributed by the transmitter and the other receiver.According to Fig. 3 we observe that by increasing the distance and time the value of γ gets closer to zero.Hence, if the distance between the receivers themselves and the transmitter is not extremely close after a certain time of observation we can assume that the barycenter is located at the center of the spherical FA receivers.In fact, authors in [21] investigated an asymptotic model for the MIMO MC system and assumed that the barycenter position could be approximated as the center of the receivers.Their justification was based on the definition of the concept and the empirical observations.From Fig. 3 we conclude that if the distance between the source is not very close, then by increasing the observation time we can assume that the barycenter is located at the center of the receiver.Studying γ and its variation contributed from the transmitter and the other receivers allow to simplify the modeling process and even skip the computation of the barycenter under certain conditions such as temporal asymptotic models [21].This analysis is one of the main contributions of this paper that allows researchers to have a tool in order to simplify their DRAFT April 27, 2023 computation while investigating diffusive MC with multiple spherical FA receivers.

B. Barycenter in SIMO
In case of p FA receivers in the channel, we can write the position of the barycenter B i of R i , as the weighted sum of contributions from the transmitter and other receivers (similar to (13)) The effect of the transmitter on the barycenter of R i is The effect of R j on R i is Accordingly the coefficients ζ (i,i) and ζ (i,j) are and

IV. SIMULATION AND RESULTS
In this section we verify the proposed model with the data obtained from PBS.In order to show the advantage of knowing the barycenter we compute the case where most of the available works assumed that the negative source is simply in the centers of the receivers.Some of the Fig. 4: A diffusive MC system with two FA receivers centered at points C 1 and C 2 .
simulation parameters are given in Tab.I, and they are borrowed from [15], except for PBS step time that we reduced even by one more order of magnitude compared to the one used in [15] to have better precision in the simulation results.We demonstrate the position of the barycenters obtained from the PBS in different scenarios and compare them with the analytical model.Moreover, we compare the expected cumulative number of absorbed molecules by the receivers with the data obtained from PBS.All the results from PBS are averaged over 100 trials and the simulation step time is 10 −7 s.The transmitter is located at the center of the coordinate.In order to simplify the visualizations we assume that all the receivers are located in xy-plane.However, the channel is 3D and the molecules move in all directions.Fig. 4 depicts the simulation scenario of SITO system.The center of R 1 denoted by C 1 is fixed on X axis at coordinate d (C 1 ,T ) , 0, 0 .
The expected cumulative number of absorbed molecules by R  The red line is drawn based on (10) while the barycenter is computed analytically as the main contribution of this paper.The green dash-dot line is when the barycenter position is substituted with the center of the receivers.We can observe an accurate match between the PBS and the  analytical results.The variation of the PBS data in case of large distances is just due to the stochasticity of the Brownian motion and considering the range of variation it can be neglected.
It can be seen that in case of the line of sight blockage of R 1 by R 2 the distribution of the particles changes severely and consequently the assumption of considering barycenter at the center of the sphere is no more valid.This observation confirms our discussion in Sec.III about the behavior of the γ when the distribution of the absorbed particles over the surface of the receiver is not uniform.The non-uniformity of the particles can be due to the closeness of a source to the receiver or the time of observation.
In Fig. 6 we depict a similar output as Fig. 5 when the radius of the receivers are not the same.In this case we assumed that the radius of the receiver R 1 is 1.2 µm and the radius of the receiver R 2 is 0.7 µm.Even in this case we observe a good match between the analytical and PBS results.The legend and markers are the same as the ones described in Fig.We observe that even in this critical case the analytical barycenters (squared markers) are very close to the empirical ones (diamond markers).
In Fig. 8a we demonstrate the variation of both analytical and empirical barycenters by The legend of the figure is the same as the one in Fig. 7.In Fig. 8a the red stars are the center of R 2 , while the red squares are the analytical barycenters, and the red diamonds are the empirical barycenters.We can observe that when Ω = 0°the attraction effect from the transmitter is strong since the distance between the receiver and transmitter is not that long.We can see that as R 2 moves behind the R 1 (the cyan star) with respect to the transmitter (the purple circle), the barycenter moves towards the center of the receiver R 2 .When Ω = 180°we see that the barycenter converges to the center of the receiver since the distance between R 2 and the transmitter increases.The repulsion effect of R 1 on the barycenter of R 2 is negligible in this case compared to that of the transmitter's effect.and the empirical barycenters (the diamond markers).We deliberately designed this configuration of receivers because we found that in the extreme cases when receivers are very close to the source or one another the estimation of the barycenter and consequently the cumulative number of absorbed molecules becomes difficult.Hence, we put R 1 very close to the source.Moreover we considered the position of R 2 such that for small angles of α it blocks the line of sight of R 4 and for the high value of α it goes to the shadowing area of R 1 and R 5 .Hence we believe that the proposed scenario is a good test to challenge different aspects of the proposed model.
In Fig. 10 we plot the cumulative number of absorbed molecules by the five receivers at t = 2 s.
The legend of the figures is the same as the one explained in Fig. 5.Each subplot corresponds to the cumulative number of absorbed particles by one of the receivers.The horizontal axis is based on the angle α in degrees.In all subfigures we can observe a perfect match between the data obtained from the empirical barycenter and PBS results.This observation once again approves the idea behind the concept of negative point source and its ideal position, which is the barycenter.
Furthermore, the red curve captures the dynamic of the empirical data and shows a very good estimation of the cumulative number of absorbed molecules.There are slight differences for receivers' observations, but given the value of the mismatch it can be considered tolerable in most practical cases.Moreover, the green dash-dot line shows a considerable deviation with respect to the PBS.This observation demonstrates the importance of barycenter analysis and consideration otherwise it is not recommended at all to neglect the barycenters and substitute them with the center of the receivers.

V. CONCLUSION
The absorption effect of fully absorbing (FA) receivers can be deduced in terms of a negative point source when there are multiple FA receivers in the channel and the best position to locate this negative point source is the barycenter.Barycenter is defined as the spatial average of the particles that are absorbed by the corresponding receiver.Localizing the barycenter is complex in diffusive MC systems because the presence of receivers brings interaction among them due to the FA characteristic.In this paper we focused on modeling the position of the barycenters in a DRAFT April 27, 2023 Fig. 10: Cumulative expected number of molecules N i , i ∈ {1, 2, 3, 4, 5} by all five receivers according to topology depicted in Fig. 9. Receiver R 2 's position varies according to the angle α, which is the relative angle between the transmitter located in the origin of the coordinate and the center of R 2 .Simulation time is t = 2 s.The legend of the figures is defined the same as the ones explained in Fig. 5.
diffusive molecular communication with multiple FA receivers.The knowledge of the barycenters allows us to solve the system of equations that describes the cumulative number of absorbed molecules.First, we derived an expression that finds the barycenter in the Single Input Single Output system.Then we impose the superposition principle and inspired by Newton's universal gravitational law we equalize the contribution of the sources (negative and positive ones) and verify the resultant model with data obtained from particle-based simulations.The simulation scenarios were chosen to be the most challenging ones, meaning that when the source and the receivers are very close but we obtained very promising results that confirm our modeling approach.We believe that the result of this paper will allow the research community to investigate the systems with multiple FA receivers.Moreover, they can also check whether it is possible to skip even the computation of the barycenter given the simulation parameters based on the

1 )
and ζ(1,2) must be designed such that represent the contribution of the two sources (i.e. the positive source as the result of the transmitter and the negative source as the result of the other receiver) on the barycenter.The positive and the negative source here do not have the same contribution to the position of the barycenter.

Fig. 1 :
Fig. 1: Center of each ring can be formulated as R cos θ where the angle θ varies from 0 to π.

Fig. 2 :
Fig. 2: Analytical γ from (24) (the red curve) versus the γ obtained from PBS (the blue curve) at t = 2 for different distances between the transmitter and the center of the receiver.

Fig. 3
Fig.3shows γ for different distances between the source and the receiver and different times of observation.We see that by increasing the time and distance, γ tends to zero.As γ becomes smaller, according to(14) the component of the barycenter that represents the effect of the source is approximately located in the center of the receiver.

DRAFT April 27, 2023 Fig. 3 :
Fig. 3: Value of analytical γ from (24) as a function of the distance between the source and the center of spherical FA receiver and time.

1
at time t = 2 s for different Ω and d (C 1 ,C 2 ) is depicted in Fig. 5.The horizontal axis of each subfigure represents the angle Ω from 0°to 180°, and the vertical axis is the cumulative number of absorbed molecules by the DRAFT April 27, 2023

Fig. 5 :
Fig. 5: Cumulative expected number of molecules N 1 absorbed by R 1 after t = 2 s, according to the scenario of Fig. 4 with d (C 1 ,T ) = 6 µm, and R 1 = R 2 = 1 µm for various positions of R 2 identified by the distance d (C 1 ,C 2 ) = {2, 4, 6, 8, 10, 12} µm and angle Ω = [0, 180]°.The red line, "Analytical", shows the cumulative expected number of absorbed molecules by R 1 based on (10) by calculating the barycenters positions analytically.The blue line, "Empirical", is using the same equation as the red line but the position of the barycenters are calculated empirically from the distribution of the particles over the receivers during the PBS.The dotted line with black squares, "PBS", is the cumulative number of absorbed molecules by R 1 .The green dashdot line, "Centered", is the value of N 1 when we assume the barycenter is located at the center of the receiver.

Fig. 8 :
Fig. 8: (a) Positions of the analytical (squares) and empirical (diamond) barycenters, B 2 , corresponding to different positions of the R 2 .The pentagram indicates the center of R 2 .Receiver R 2 revolves around R 1 at distance d C 1 C 2 = 4.(b) Magnified position of the empirical (diamond) and estimated (square) barycenters, B 1 .

Fig. 8b shows
Fig. 8b shows the variation of the barycenter B 1 corresponding to different positions of R 2 shown in Fig. 8a.As the position of the R 2 changes clockwise the position of the barycenter B 1 also changes clockwise.Note that Fig. 8b is a zoomed version of a small area inside R 1 and the difference between the analytical barycenters (the blue squares) and the empirical ones

TABLE I :
Simulation parameters