The Method of Fictitious Negative Sources to Model Diffusive Channels With Absorbing Boundaries

This paper presents an approach to address the diffusion equation in scenarios involving multiple absorbing boundary conditions, commonly found in diffusive molecular communication (MC) channels. Instead of using multiple mirror images of the source, fictitious sources with time-varying release rates are introduced to replace the boundaries. This transformation enables the calculation of the expected cumulative number of absorbed particles (CNAP) by multiple absorbing boundaries with finite volume. To compute the expected CNAP, the concept of barycenter, which represents the spatial mean of particles the receiver absorbs is introduced. Substituting absorbing objects with their barycenters leads to model the CNAP in scenarios with convex geometry of absorbers. In a one-dimensional (1D) space, the proposed approach yields the same expression as the method of images for describing the expected CNAP by an absorber. However, in three-dimensional (3D) space, where using the method of images is challenging or even impossible, the proposed approach enables substituting the objects with fictitious sources and compute the expected CNAP. In 1D, an extension of this approach to the case in which one boundary exhibits an absorption characteristic while the other has zero-flux characteristic is demonstrated. This research direction is valuable for modeling channels where not all objects are particle receptors.

The Method of Fictitious Negative Sources to Model Diffusive Channels With Absorbing Boundaries Fardad Vakilipoor , Abdulhamid N. M. Ansari , Luca Barletta , Member, IEEE, Gian Guido Gentili , Member, IEEE, and Maurizio Magarini , Member, IEEE Abstract-This paper presents an approach to address the diffusion equation in scenarios involving multiple absorbing boundary conditions, commonly found in diffusive molecular communication (MC) channels.Instead of using multiple mirror images of the source, fictitious sources with time-varying release rates are introduced to replace the boundaries.This transformation enables the calculation of the expected cumulative number of absorbed particles (CNAP) by multiple absorbing boundaries with finite volume.To compute the expected CNAP, the concept of barycenter, which represents the spatial mean of particles the receiver absorbs is introduced.Substituting absorbing objects with their barycenters leads to model the CNAP in scenarios with convex geometry of absorbers.In a one-dimensional (1D) space, the proposed approach yields the same expression as the method of images for describing the expected CNAP by an absorber.However, in three-dimensional (3D) space, where using the method of images is challenging or even impossible, the proposed approach enables substituting the objects with fictitious sources and compute the expected CNAP.In 1D, an extension of this approach to the case in which one boundary exhibits an absorption characteristic while the other has zeroflux characteristic is demonstrated.This research direction is valuable for modeling channels where not all objects are particle receptors.
Index Terms-Molecular communication, diffusion, channel modeling.

I. INTRODUCTION
C OMMUNICATION and information exchange have always played an essential role since the beginning of the universe.Matter, energy, and information are the essential three elements that construct the universe.Communication extends beyond the realm of human interaction, and it encompasses the information exchange in nature.The propagation of waves has always been the core methodology for exchanging information.Electromagnetic and acoustic waves have been extensively studied, whereas particle-based propagation has received comparatively less attention in research pertaining to communication studies.One of the primary forms of propagation of particles is diffusion.In this work, we dedicate our attention to the diffusive propagation of particles and attempt to investigate the expected cumulative number of absorbed particles (CNAP) by the boundaries with absorption characteristics.Diffusion is a fundamental mechanism that enables molecule movement and interaction in biological environments.Hence, it has received the most attention among the propagation mechanisms existing in information exchange at the biological level utilizing particles.
Recent advancements in synthetic biology and bionanotechnology have catalyzed transformative prospects across diverse domains, including medicine, tissue and materials engineering, and environmental monitoring [1].Central to these strides is the burgeoning potential of cooperative interactions among bionanomachines, an achievable feat necessitating network organization and viable internode communication [2].Bionanomachines can exchange information in fluidic environments using the molecular communication (MC) paradigm, which involves the exchange of signaling particles [3].
Examples of micro-scale communication [4] include communication among cells [5] or bacteria [6], [7], while macro-scale [8] examples include the olfactory system [9] pheromonal communication among animals [10].The biocompatibility and energy efficiency of MC make it a promising information transmission paradigm in the sub-millimeter range.MC applications have been mainly focused on the health sector, where it has shown potential for smart drug delivery [11], [12], [13], health monitoring [14], and disease diagnosis [15], among other areas.By bridging the gap between theory and reality, MC holds promise for a variety of practical applications.
The presence of multiple receivers is the crucial element in the MC domain.Due to the existing interactions between the receivers and particles, which are mainly interpreted as an absorption effect, the receivers exhibit interaction among themselves.The interaction reveals itself by comparing the CNAP associated with the receiver in the presence and absence of the second receiver in its proximity.As always, to investigate and characterize the communication aspect of such a system, the knowledge of the channel impulse response (CIR) would be necessary and helpful.Another motivation in investigating the presence of multiple receivers is that such systems are closer to reality since an individual receiver alone cannot afford the complexity of the task.In natural MCs 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.

A. Related Literature
There has been interest in MC literature to move toward crowded environments where there are multiple objects with absorbing features in the physical channel.From the differential equations perspective, these objects are equivalent to the boundary conditions that must be considered mathematically.
Numerous channel models for MC with diffusive propagation of the particles considered the impact of neighboring multiple absorbing receivers [16], [17], [18], [19].In order to count differences in the number of particles absorbed in the case of two fully absorbing (FA) receivers compared to one, a preliminary investigation is conducted in [16].Authors in [17] derived the identical analytical closed-form expression for the fraction of particles absorbed over time by each receiver using the same one-dimensional (1D) environment.Of course, the point transmitter is always positioned between the two absorbers, which is a feature of the 1D space.Consideration of the two FA receivers in a three-dimensional (3D) environment is suggested in [18], from which the fraction of absorbed particles by each receiver is computed.
In an effort to deduce an analytical formulation adjusted to the complex 3D geometry with two FA receivers, an endeavor is launched in [20].This formulation tackles encapsulating the influence of a second receiver by conducting a scaling coefficient that influences the number of particles absorbed by the intended receiver up to a certain time instant.The scaling factor is determined by curve-fitting empirical data and is a function of the geometrical characteristics with certain adjustment coefficients unique to each receiver's position.The main technical drawback of the approach is that it makes the assumption that the effect of the interfering receiver can be described by a straightforward decrease in the number of particles released by the source, but in reality, the shape of the CIR can be disrupted.Additionally, it is difficult to generalize the procedure to any quantity of receivers.Another similar approach of expressing the interaction between receivers by finding constant coefficients from the empirical data can be found in [21].The main focus of [21] was investigating a modulation technique, but at first, they looked for a model to capture the effect of multiple receivers in the channel.Authors in [22] followed a similar approach as in [20].They modeled the CIR using a function similar to the single-input single-output (SISO) system response by applying curve-fitting algorithms.Specifically, they used control coefficients over the SISO model to describe the MIMO system.The control coefficients were selected to comprehend system characteristics and be time-independent.Then, nonlinear regression models were used to fit simulation data.
Authors in [23] derived 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 transmitter, as well as between the receivers themselves, are sufficiently large.Recently, a semi-analytical model proposed in [24] for arbitrary geometry of channel included a first order chemical reaction.By employing the numerical method of moments, they transform the problem of concentration Green's function (CGF) derivation from the CGF linear integral equation (CLIE) into an inverse matrix problem.
Previously in [25], we 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 receivers in the physical channel was taken into account by describing them as point sources of negative particles.However, replacing a volumetric object with a point requires extra attention about where to locate the point.We proposed the concept of barycenter as the spatial mean of the particles that are absorbed by each receiver, and developed an empirical expression to describe the position of the barycenter at the time of observation t = 2 s by applying curve fitting and finding an empirical expression.We 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 surface-point of the spherical FA receiver.The surface-point is defined as the closest point on the surface of the receiver to the transmitter.We proposed a parameter γ that varies between 1 and 0. When γ is equal to 1, it means that the spatial average of the absorbed particles are all concentrated at the surface-point of the receiver, while γ equal to 0 means that the absorbed particles are distributed around the receiver, and their spatial average coincides with the center of the receiver.After [25], authors in [26] proposed a channel model similar to the concept proposed previously based on the concept of missing particles by the absorbing boundaries, and they obtained solutions with good accuracy when there is no strong shadowing effect in the system.
Recently, we focused on modeling the barycenter in 3D space with multiple FA receivers in [27].We derived a heuristic analytical expression that serves to determine the location of the barycenter within a spherical FA receiver as a function of time.Furthermore, the resulting model affords the opportunity to gain insights into the specific conditions under which simplifying assumptions can be employed.In certain cases, it becomes feasible to circumvent the intricacies of barycenter computations and instead adopt an approximation where the barycenter is assumed to be positioned at the center of the receivers, which means that the parameter γ from [25] is set to zero.To demonstrate the importance of computing the barycenter, we also considered the case where the negative source is assumed to be at the center of the receivers.Comparing the expected CNAP from two cases with the results from the particle-based simulation (PBS) depicted the importance of computing the barycenter.

B. Motivation and Contribution
Most of the papers discussed in Section I-A attempt to describe the system model from the probabilistic approach.However, in this work, we demonstrate and analyze a method to solve the diffusion equation with different boundary conditions (absorbing and zero-flux conditions).In certain scenarios, we demonstrate that our proposed method yields identical results to the method of images [28] for solving the diffusion partial differential equation (PDE).Moreover, for scenarios where the method of images fails, we show how the proposed approach gives accurate approximation to the solution of the diffusion PDE.First, we discuss the method of images.We demonstrate that this method requires considering multiple mirror images of the source with respect to the boundary.The absorbing feature of the boundary results in the reduction of the particles from the environment.This observation inspired us to investigate the possibility of replacing the absorbing boundaries with fictitious sources that release negative particles.
The advantage of using fictitious sources in place of the method of images is twofold: • Introduction of a finite number of fictitious sources, equal to the number of boundaries, compared to the introduction of infinite replicas of the main source as suggested by the method of images; • In the presence of 3D convex bodies, the position of the fictitious source is always inside the volume of the body, thus enabling us to formulate a model that approximates the expected CNAP in the receivers because we are confident that the barycenter is inside the receiver (i.e., equivalent to a finite space/domain).Finally, for scenarios where there is also a zero-flux boundary condition in addition to the absorbing boundary, we obtain the expected CNAP in 1D space.We demonstrate that our proposed approach results in the same solution as that obtained from the method of images.This last example demonstrates the power of the proposed method to open the doors towards modeling complex environments when not necessarily all the boundaries are the receptors of the particles.
We first compute the CNAP from the PBS simulation.Then, we compute the barycenter according to its definition from the distribution of the absorbed particles over the surface of the receivers and substitute the coordinates of the barycenters into the proposed model as the locations of the fictitious negative sources.By comparing the CNAP from PBS and the proposed method, we verify the proposed system model and the concept of barycenter.We deliberately considered a complex topology of the receivers' configuration where five spherical receivers with different sizes are close to one another and to the source.To the best of our knowledge, such a complex scenario has never been considered in the MC literature to verify the channel models.One should note that the arrangement of the spherical boundaries in a crowded biological environment is unpredictable; hence, to ensure the generalization of the concept, we seek to asses a complex scenario with a dominant shadowing effect.Moreover, we devise another simulation scenario to demonstrate that there are cases where we can skip the computation of the barycenter and replace it with the center of the receivers.The latter case, can be regarded as the simple scenario.We would like to mention that compared to our previous related works [25], [27], in this paper we focus on using the technique of fictitious time-varying negative source instead of focusing on modeling the barycenter point.To be more specific, we start with classical method of images and provide a novel way of looking at the system of equations describing the interaction among the boundaries.To the best of our knowledge, none of these demonstrations and discussion have been pointed out previously in the literature.

C. Outline
The rest of the paper is organized as follows.Section II discusses the method of images and then obtains the CNAP in 1D space, starting with a single absorber and extending to the two absorbers.In the same section, we also introduce the proposed modeling approach of replacing the absorbers with negative sources and demonstrate that it results in the same solution as the one attained from the method of images.In Section III, we take advantage of the fictitious source method and demonstrate how to model the system when there are multiple spherical receivers in a channel with 3D geometry.Section IV depicts the power of the proposed approach in taking into account the effect of the zero-flux boundary.It is proved that in 1D space when there is a zero-flux boundary and an absorbing boundary, the CNAP from the method of images is the same as the one obtained from the proposed method.In Section V, we depict and discuss the simulation results.First, we show that the idea of the barycenter is valid and that it is the best position where to place the negative source.Then, we place receivers randomly in space with a certain minimum distance from one another and show that even approximating the barycenter as the center of the receivers results in a promising accuracy in such a scenario.Finally, Section VI provides the conclusion of this paper.

II. 1D SPACE WITH ABSORBING BOUNDARIES
In this section, we discuss the method of images first to obtain the solution of the diffusion equation and then to attain the rate of the absorbed particles by the absorbing boundary.Then, we describe the other possible interpretations to represent the model and show that we can achieve the same solution in terms of the expected rate of the absorbed particles by the boundary as the one attained from the method of images.

A. 1D Model With Single Absorber
Before initiating the investigation on modeling the system that contains multiple boundaries, let us start with the case of a single boundary that allows us to introduce the method of images briefly.We start the modeling of the 1D diffusion in space with an absorbing boundary condition located at the origin of the coordinate system (x = 0) and a dimensionless source with an impulsive release feature at x = L. Diffusion equation is a parabolic PDE that connects the first-order Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.temporal variation to the second-order spatial variation of the desired variable with a diffusion coefficient D. Here, we consider the normalized concentration C(x, t) that is the number of particles per unit of time and space divided by the total number of released particles by the source.The normalized concentration C obeys the diffusion equation with boundary and initial conditions where δ(•) is the Dirac delta function, IC and BC stand for the initial condition and boundary condition, respectively.One of the common approaches to devise a solution for the diffusion equation is the method of images.In the method of images, we can separate the effect of the source and boundaries.The solution of the diffusion equation excited by an impulsive source in an unbounded environment is described by a Gaussian function.Then, to satisfy the presence of an absorbing boundary condition, an image of the source with the opposite effect at x = −L needs to be considered to force the concentration to be zero at the boundary, see Fig. 1.We can write the solution with an absorber as the superposition of the two sources with opposite signs at two specular locations with respect to the absorber.Hence, the solution of (1) subject to (2) can be written as where x ≥ 0 and t > 0.
The fraction of free particles at time t (fraction of particles not absorbed by the time t with respect to the total number of released particles) in the space can be found by the spatial integration of the normalized concentration in the free space domain [0, ∞] (see Appendix A) In this specific example, if we multiply 1 − F(t) by the total number of particles instantaneously released from the source, we obtain the expected CNAP by the absorber.The rate of change in normalized concentration of the particles in free space can be obtained by taking the derivative of ( 4) with respect to time.Obviously, adding a negative sign is required to allow us to achieve the rate of fraction of particles being absorbed from the absorber's perspective.We can write the absorption's rate of normalized concentration at the boundary with distance L from the source in 1D space as If we assume the particles' movements are independent and the impulsive source is able to release N T particles, then we define where n(t) is the rate of absorption by the corresponding boundary.

B. 1D Model With Two Absorbers
After defining the rate of absorption in the case of a single absorber in the channel, we want to obtain the rate of absorption when two absorbers are present.We demonstrate the solution from the method of images and then show that the approach proposed in this work results in the exact same expression.Similar to the previous example with one absorber, from the method of images, the goal is to force the concentration to be zero at the absorbing boundaries.With reference to the geometry shown in Fig. 2, to neutralize the concentration at x = L we place a negative source at x = 2L.Note that the tail of the Gaussian kernel of this negative source also affects the concentration at boundary x = −L.Then, to compensate this side effect, we are forced to place a positive source at x = −4L, i.e., in the symmetrical viewpoint of the negative source at x = 2L with respect to absorber at x = − L. With similar reasoning, the same methodology must be applied to satisfy the presence of the other absorber.Hence, the concentration between the two absorbing boundaries can be seen as a superposition of the effect of infinite negative and positive sources.These sources are essentially mirror images of the original source with respect to the boundaries, as depicted in Fig. 2.
The solution of the diffusion equation, according to the method of images, can be written as where x ∈ [−L, L] and t > 0. The rate at which the normalized concentration of absorbed particles changes is equal in magnitude but opposite in sign to the rate of reduction in the normalized concentration of free particles.The reduction rate can be obtained by taking the spatial derivative of ( 7) at x = L (we are interested in the absorption rate at the boundary located in x = L) multiplied by D according to which is based on (1).Hence, we write the rate of absorption at the boundary x = L as After obtaining the result with the interpretation of mirror images, we would like to discuss and introduce another interpretation to obtain the absorption rate of the boundary in the same physical environment.Authors in [18] provided a beautiful proof for a system of equations to approximate the interaction between the two absorbing boundaries.They partitioned the probability of the events a particle can take.Inspired by their resultant equation we provide a novel viewpoint of looking at the integrated system of equations to capture the interaction among the boundaries in general.Instead of considering the image reflection of the main source, the other way of looking at the problem is to replace the absorbing boundary with sources that release negative particles.Hence, one can describe the model in terms of a coupled system of equations as where is the convolution operation.Each equation in (10) describes the rate of absorbed particles at the specific boundary.The rate of absorption at boundary x = L and x= − L is indicated by n 1 (t) and n 2 (t), respectively.The rate of absorbed particles at x = L can be seen as the contribution of the real (positive) source at distance L from the absorber and the fictitious (negative) source to capture the effect of the other absorber at x = − L. The contribution of the positive source is the impulse response of the absorption rate as if there were not the second absorber around.Due to the absorption feature of the left side boundary that is equivalent to reducing the particles from the environment from the perspective of the right side receiver, a fictitious source is assumed that releases negative particles.The release rate of negative particles by the fictitious source is exactly equal to the rate of absorbed particles by the absorber.It is important to note that the system of equations in ( 10) is consistent with [18, eq. ( 29)], which explains a similar system in its integral form that was derived using a probabilistic approach.In Appendix B, we show that the solution of ( 10) is equal to (9).After ensuring that our system of equations allows us to capture the interaction between the absorbers in 1D, we would like to step further for complex boundaries in higher dimensions.

III. 3D SPACE WITH SPHERICAL ABSORBING BOUNDARIES
In reality, objects have a finite volume that does not dissect the space, and the 3D geometry of the space must be considered.We desire to find an approach that allows us to obtain the rate of absorbed particles and, consequently, the expected CNAP associated with the absorbers when multiple objects with absorbing surfaces are present in the environment.These objects can be recalled as receivers from the communication perspective.Let us assume a 3D space with spherical absorbers.This geometry is consistent with real bacteria with coccus cell morphology (e.g., Staphylococcus species, Streptococcus species, etc.), i.e., round-shaped bacteria whose diameter vary in the range [0.25, 1] μm [29].
In this case, we cannot rely on the image method to create multiple mirror images of the source with respect to the absorber, because the absorbing boundary condition no longer dissects the space.Intuitively, the dissection of the space by the boundary allowed us to have mirror images of the source, but in the case of finite volume objects, devising the image of the source with respect to the boundary is complicated and unimaginable.However, by following the proposed method (10), we can substitute the receivers with fictitious negative sources of particles and solve the system of equations.According to [30], the rate of absorbed particles by the spherical absorber for t > 0 when there is no other absorber in the space reads where d (C,T ) is the distance between the center of the spherical absorber C and the source/transmitter T .The radius of the receiver is denoted as R.
In the case of 1D space, we simply replace the absorber with the negative source.This substitution is consistent with the definition of the barycenter.The barycenter is defined as the spatial mean of the particles that hit the surface of the absorber.Obviously, in 1D space the mean of the position of the particles that hit a certain absorber is exactly where the absorber is located.Unlike the 1D space, in 3D with finite volume boundaries, replacing the volumetric objects (receivers) with negative sources is challenging.
In [25], the concept of the barycenter was investigated to solve the system of equations and obtain the expected CNAP by the spherical receivers.We showed that to obtain the expected CNAP for a fixed time t = 2 s, the best position to place the fictitious negative source is the barycenter of the absorber.The difficulty in solving the system of equations in 3D space with finite volume receivers is that the barycenter of the receivers varies with time because the distribution of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the particles that are hitting the surface of the receivers varies with time [27].
Even though the position of the barycenter is expected to vary with time, we separate the dynamic in the variation of the position of the barycenter and the expected number of absorbed particles.To model the barycenter exclusively, one can rely on the proposed expressions in [25], [27] or implicitly compute the distribution of the absorbed particles over the surface of the boundary and compute the mean according to the barycenter's definition.Demonstrations regarding the barycenter dynamic can be found in the aforementioned references.Note that the barycenter is the best position to place the negative source to obtain the expected CNAP.Hence, while applying the Laplace transform to the system of equations to solve it, we treat the position of the barycenter as an unknown parameter independent of time.We separately focus on where to put the barycenter and model it.In Section V, we demonstrate that if we compute the barycenter at every time instant from the PBS data, which is the ground truth of the position of the barycenter, and substitute it in the integral form of the system of equations or its closed-form solution, the result accurately describes the expected CNAP by the receivers.Hence, solving the system of equations and modeling the position of the barycenter as two separate problems is a reliable solution.Decoupling the two problems of solving the system of equations and of placing the position of the barycenter is another advantage that the fictitious source method provides us.
We start writing the system of equations to express the rate of absorbed particles by the receivers with the final goal of computing the CNAP.In the integral form, the system of equations is where B 1 and B 2 are referring to the barycenter points of the receivers R 1 and R 2 , d (C 1 ,B 2 ) is the distance between the center of spherical receiver R 1 and barycenter B 2 , and d (C 2 ,B 1 ) is the distance between the center of spherical receiver R 2 and barycenter B 1 .However, we take the derivative (12) and write it in terms of interaction between the rate of absorption by the receivers.In (12), N 1 (t) and N 2 (t) are the expected CNPAs corresponding to the two receivers.
The first receiver sees the other one as a negative point source located at its barycenter.The absorption effect of the receivers from one another perspective is seen as a reduction of the particles from space.This feature is replaced with a negative source such that its rate of releasing negative particles is equal to the rate of absorbing particles.Finally, we take the derivative of ( 12) with respect to time and neglect the dependency of the barycenter on time for the moment.The reason behind this assumption has roots in the time scale separation of the two variables (i.e., rate of the absorption and the barycenter).According to the definition of the barycenter, it is the temporal integration of the distribution of the particles that are being absorbed by the boundary from time zero, while the rate of the absorption considers the rate of particles being absorbed for the moment.Hence, these two concepts are not from the same temporal scale.For this reason, we assume that the Barycenter position can be seen as time-invariant within short time intervals with respect to the variation of the absorption rate.However, this should not imply that the barycenter is time-invariant in general because, as we showed in [27] it depends on the particles' distribution through time.Accordingly, this strong assumption allows us to decouple the problem of locating the barycenter from the problem receivers' interaction with each other to some extent.The expected rate of the absorbed particles by the spherical receivers is described by , (13) The solution of ( 13) is ( 14), shown at the bottom of the page (see Appendix C).
In general, we are interested in the CNAP.This concept is equivalent to taking the integral of the rate of absorbed particles in the time domain starting from the time zero.The integral operation is equivalent to divide the rate by s in the Laplace domain.
The CNAP by the spherical absorber is ( 16), shown at the bottom of the next page.The introduced approach can Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
be extended to multiple receivers in the environment by expanding the superposition principle and the coupling effect of the boundaries on each other (please check [27] for the generalization of the model); however, the closed-form solution to the system of equations is not obtained, and we rely on numerical integration.In the case of convex-shaped receivers, the barycenter is always inside the volume of the receiver based on its definition.Hence, it is possible to investigate and apply the assumptions to neglect the computation of the barycenter as hypothesized in [27] and [31].The method of fictitious sources allows us, instead of looking at the whole space to place the mirror images of the source with respect to the boundaries, to deal with a finite space inside the volume of the absorber to place the negative source.This feature allows us to find the asymptotic position of the negative source subject to the geometry of the receivers and the topology of the channel.

IV. 1D SPACE WITH ABSORBING AND ZERO-FLUX BOUNDARIES
In general, different boundaries can exist in the biological channel where some can be interpreted as the receiver or their presence affects the intended receiver's observation.The most commonly studied one is the transparent receiver, which is called the perfect monitoring receiver in biophysics [32].In transparent receivers, the boundary does not interact with the particles; hence, it is just a volumetric integration in time over the concentration of particles in space.An interesting problem that requires more attention in MC studies, to the best of our knowledge, is when some objects in the environment are not the receptors of the particles.This assumption is another step to allow researchers to investigate complex and crowded environments.The absorption feature of the receivers actually comes from the assumption that the reaction rate between the receiver's surface and the particles goes to infinite.So, the preliminary assumption is that the receiver is capable of having the desired reaction with particles, but this assumption does not always hold in reality.There are surfaces that do not react with the particles.A common assumption to model the behavior of the particles at this type of boundary is the zero-flux characterization.The zero-flux boundary indicates that the variation of the concentration with respect to time is zero where the boundary is located.Hence, the particles cannot pass the barrier or be absorbed.The zero-flux boundary can be described as Let the absorbing boundary be located at x = L and the zero-flux boundary be at x = − L.
From the method of images, in order to satisfy the effect of the zero-flux boundary at x = − L we place a positive source at x = − 2L to create the particles accumulation effect at the boundary x = − L. Subsequently, we place a negative source at x = 4L to neutralize the contribution of the aforementioned imaginary positive source at the absorbing boundary in x = L.Then, to describe the effect of the absorber, we place a negative source at x = 2L, and to capture the effect of this negative source at the zero-flux boundary, we need to add a negative source at x = − 4L.Then, to neutralize the effect of the negative source at the absorbing boundary, we place a positive source at x = 6L and so on (see Fig. 3).In the end, we can describe the normalized concentration of the particles in terms of an infinite series as According to the definition and discussion from (8), we obtain the rate of absorbed particles by the absorbing boundary as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
This problem can also be formulated from the coupled system of equations approach by replacing the boundaries with fictitious sources.Actually, the zero-flux boundary can be seen as a positive source from the perspective of the absorber at the distance 2L from it, because the zero-flux boundary is performed like a solid wall that a particle cannot cross and remains in the free space i.e., x ∈ [−L, L].This behavior is equivalent to assuming that a particle has passed the zero-flux boundary, and the fictitious source generates a new particle.On the other hand, the absorber can be seen as a negative source from the perspective of the zero-flux boundary due to its feature in reducing the particles from space.
First, from the perspective of the absorber, there is a source at the center of the coordinate system at distance L from it.Then, the zero-flux boundary effect is described as a source that releases positive particles with a time-varying rate.A similar intuition holds to describe the rate of particles passing the absorbing boundary.In this case, the system of equations requires a slight change compared to (10) as The solution of ( 20) is ( 19) (see Appendix D).We showed that the new approach allows us to simply obtain the expected rate of the absorbed particles and, consequently, the expected CNAP by the receiver in 1D space.Consideration of the zero-flux boundary is another part of investigating diffusive channels, as not all the boundaries are receptors of the particles.Similar to the previous example, the boundaries here are dissecting the space, and the positions of the fictitious sources are fixed.However, we found that in 3D space with finite volume objects featuring zero-flux properties, the concept of barycenter does not results in an accurate solution.Hence, a new concept is required to be defined for future works to allow us to investigate complex environments containing finite volume objects with both absorbing and zeroflux surfaces.

V. SIMULATION AND RESULTS
In this section, we verify the concept of barycenter and the accuracy of the proposed approach for obtaining the CNAP by the spherical receivers over time when the barycenter is known.
For the simulation scenario, we consider the most complex receiver topology compared to those found in the existing literature on channel modeling with multiple FA spherical receivers.This is done to demonstrate the validity of the barycenter concept and the effectiveness of our proposed approach.The propagation medium has a diffusion coefficient D of 79.4 μm 2 /s, which is equivalent to a water-based diffusive environment in the human body [33].Note that D affects the propagation speed of the particles.Hence, the particles' distribution on the boundary's surface rapidly converges to a steady state, which implies that the barycenter reaches a fixed value according to its definition.The ground truth of our theoretical evaluation was obtained from the PBS data.The simulation step time is 10 −7 s.Every simulation scenario is repeated for 100 times, and then we take the average of the CNAP by each receiver.To verify our system model, we compute the barycenter of each receiver according to its definition from the PBS data, and then substitute it in the extended version of (12).Due to the simplicity of extending (12) to more than two receivers and in order to save space, we skip writing the extended form, although it is straightforward from intuition.It is important to note that no closed-form solution has been devised yet for the extended version of ( 12).However, it can be solved simply by using a numerical integration method.
Figure 4 depicts the first simulation scenario we considered in this paper.The source is located at the center of the 3D coordinate system.Receiver R 1 is at (3 × 10 −6 , 0, 0) with radius 1 μm.Receiver R 2 is at (7 × 10 −6 , 0, 0) with radius 2 μm.Receiver R 3 is at (0, 3 × 10 −6 , 0) with radius 2 μm.Receiver R 4 is at (0, 0, −5 × 10 −6 ) with radius 3 μm.Receiver R 5 is at (0, 7 × 10 −6 , 0) with radius 1 μm.We deliberately considered a case where receivers are in the shadowing area of the other receivers.Moreover, in one case, we put the smaller receiver in the shadow area, while in the other case, we located the bigger receiver in the shadow area.Furthermore, we also located a big receiver close to the source and other receivers.Despite previous works, we considered the case where receivers are very close to the source to verify the concept's validity.Note that since the receivers are very close to the source, the CNAP corresponding to each receiver reaches its asymptotic value in a short time.In order to better demonstrate the accuracy of the model, we consider a simulation time of 0.6 s.Please note that, barycenter modeling is a complicated task that needs proper attention.For instance, we would like to draw the reader's attention to compare [26, Fig. 5(a)] and [27, Supplemental Items, Fig. 47].Both figures demonstrate the CNAPs of the absorbing receivers corresponding to the same scenario, while the difference is notable.This difference between the obtained results goes back to different analytical approaches employed by the authors.
Figure 5 shows the expected CNAP by the receivers.The solid lines represent the results obtained from the PBS.The Cumulative number of the absorbed particles by the receivers according to the simulation scenario shown in Fig. 4.
color of the curves is chosen to be the same as the spheres in Fig. 4. The circular markers correspond to the results obtained by solving the system of equations and substituting the barycenter computed from the PBS data.We can see that there is an excellent match between the results in such a complex scenario.This numerical experiment validates the idea of the barycenter to express the position of the negative source.Note that without the concept of the fictitious negative source, we could not have reached this stage to model the receivers' observations.At this stage, researchers can focus on modeling the barycenter or find scenarios where they can approximate its position.By relying on the method of images, one could not have reached this stage to attack the problem of multiple objects in 3D space.Moreover, to show the importance of the role of the barycenter in modeling, we substitute the barycenter with the center of the receivers to show the difference (see dotted lines with squared markers).We observe that by misplacing the negative sources in the center of the receivers, the receivers' observation becomes negative, which does not have any physical interpretation.In the realm of fictitious sources, this negative observation can be interpreted as the negative particles absorbed by the fifth receiver being more than the positive particles.
Figure 7 reports the CNAP by the receivers according to the simulation scenario shown in Fig. 6.The solid lines correspond to the CNAP obtained from the PBS.The dotted lines with squared markers correspond to the solution of the  Cumulative number of the absorbed particles by the receivers according to the simulation scenario shown in Fig. 6.
system of equations if we neglect the role of the barycenter and simply replace it with the center of the receivers.We can observe that in this scenario, disregarding the barycenter is a valid approximation because the receivers are not extremely close to one another or to the source.

VI. CONCLUSION
This paper introduces and discusses a new approach to solve the diffusion equation with absorbing and zero-flux boundary conditions.Inspired by the method of images to solve the diffusion equation, instead of considering multiple mirror images of the source with respect to boundaries, we replace the boundaries with fictitious sources featuring a timevarying release rate.
We showed that in 1D space, the proposed approach results in the same expression as we can obtain from the method of images when two absorbers are present in the channel.Then, we discussed a similar problem in 3D geometry while the absorbers do not dissect the space and they have a finite volume.In this case, the position of the corresponding fictitious source varies with time because the spatial mean Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
of the particles hitting the surface of the absorbers is timevarying.However, at this point, the advantage of the proposed approach comes in handy.In the method of images, we consider reflections of the source with respect to the boundaries, and due to the interaction between the boundaries, this results in considering infinite mirror images of the source, but in the proposed approach, we replace the boundaries with a single fictitious source.Hence, in an environment with a finite number of absorbers, we have a finite number of fictitious sources.
In computing the expected CNAP, the best place to locate the fictitious source is the barycenter, which is always inside the volume of the absorber with convex geometry.Knowing that the position of the source is inside the volume of the absorber allows us to investigate scenarios where we can approximate the position as the center of the spherical object, especially when we are interested in the CNAP after a long time passes from the release moment.
Lastly, we discussed a scenario that includes a zero-flux boundary and an absorbing boundary.We demonstrated that for 1D problems, the proposed interpretation results in the exact same solution that can be obtained from the method of images.However, extension of the fictitious source method to 3D space has not been solved yet, because the position of the fictitious positive source that is supposed to be replaced with the zero-flux boundary is unknown and worth investigating in the future.In general, we believe that this new approach allows researchers to obtain the expected CNAP in a physical channel where the particles' movement is governed by the parabolic PDE (diffusion equation) and multiple objects are present in the channel.Furthermore, it is worth investigating the applicability of the proposed approach in solving hyperbolic PDEs (wave equation).

APPENDIX A
The fraction of the particles in free space by the time t with respect to the total particles released by the source instantaneously is the spatial integral of the concentration over the free space domain Let us change the variable then where . We can also write the fraction of the absorbed particles with respect to the total number of released particles as 1 − F (t), which is erfc( L √ 4Dt ).

APPENDIX B
Let us take the Laplace transform of (10) to get rid of the convolution operation and convert it into multiplication.Hence, we can write (10) in Laplace domain as The solution of this system of equations for n 1 (s) is We substitute 1 1+f (2L,s) with its power expansion series The Laplace transform of f (L, t) is Thus n 1 (s) is The inverse Laplace transform of n 1 (s) is (−1) j (2j + 1)L √ 4Dπt 3 e − ((2j +1)L) 2 4Dt Let us decompose the series in terms of even and odd increments of the summation.
We multiply and divide (30)  .(31) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
We change the sign of the increment in two summation series .(33) In the end, ( 9) and ( 33) are exactly the same.The solution of the system of equations (34) is In a similar way as before, we write the power expansion series equivalent to the inverse of the denominator in fraction (35) ( The time domain solution of the rate of absorbed particles by the receiver one is (14).

APPENDIX D
We follow the same procedure as of Appendix B. We apply the Laplace transform to (20) The solution of this system of equations for n 1 (s) is Substituting 1 1+f 2 (2L,s) by its equivalent power expansion series (41) (−1) j f (L, s)f 2j (2L, s) We use the Laplace transform ( 27) (44) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 1 .
Fig. 1.Image method representation of the solution to the diffusion equation with a single absorber placed at x = 0 (the vertical black dashed line) in 1D.

Fig. 2 .
Fig. 2. Image method representation of the solution to diffusion equation with two absorbers (vertical dashed black lines) at x = L and x = −L in 1D space.

Fig. 3 .
Fig. 3. Image method representation of the solution to diffusion equation with an absorber (vertical dashed black line) at x = L and a zero-flux boundary (vertical solid black line) at x = − L in 1D space.

Fig. 5 .
Fig. 5.Cumulative number of the absorbed particles by the receivers according to the simulation scenario shown in Fig.4.

Fig. 6 .
Fig. 6.Topology of the simulation scenario that is possible to skip the barycenter computation.

Fig. 7 .
Fig. 7.Cumulative number of the absorbed particles by the receivers according to the simulation scenario shown in Fig.6.