Microfluidic QCSK Transmitter and Receiver Design for Molecular Communication
11 Microfluidic QCSK Transmitter and ReceiverDesign for Molecular Communication
Dadi Bi,
Student Member, IEEE,
Yansha Deng,
Member, IEEE,
Abstract
The design of components with molecular communication (MC) functionalities can bring an oppor-tunity to enable some emerging applications in fields from personal healthcare to modern industry. In thispaper, we propose the designs of the microfluidic transmitter and receiver with quadruple concentrationshift keying (QCSK) modulation and demodulation functionalities. To do so, we first present an ANDgate design, and then apply it to the QCSK transmitter and receiver design. The QCSK transmitter iscapable of modulating two input signals to four different concentration levels, and the QCSK receivercan demodulate a received signal to two outputs. More importantly, we also establish a mathematicalframework to theoretically characterize our proposed microfluidic circuits. Based on this, we first derivethe output concentration distribution of our proposed AND gate design, and provide the insight intothe selection of design parameters to ensure an exhibition of desired behaviour. We further derive theoutput concentration distributions of the QCSK transmitter and receiver. Simulation results obtained inCOMSOL Multiphysics not only show the desired behaviour of all the proposed microfluidic circuits,but also demonstrate the accuracy of the proposed mathematical framework.
Index Terms
Molecular communication, microfluidics, signal processing, chemical reactions, AND gate, QCSKmodulation and demodulation.
I. I
NTRODUCTION
Over the past few years, molecular communication (MC) has attracted increasing atten-tion as it can wave revolutionary and interdisciplinary applications ranging from healthcare,to industry, and military [1], [2]. This inspired a bulk of research centering around theoreticalcharacterizations of MC, such as transmission schemes [3], [4], propagation characterizations[5], [6], reception models [7], [8], and detection strategies [9], [10]. To ensure successfulinformation transmission, signal processing units are envisioned to be essential componentsfor MC transmitter and receiver to facilitate modulation/demodulation and coding/decoding
D. Bi and Y. Deng are with the Department of Engineering, Kings College London, London, WC2R 2LS, U.K. (e-mail: { dadi.bi,yansha.deng } @kcl.ac.uk). (Corresponding author: Yansha Deng). a r X i v : . [ c s . ET ] M a y functionalities. However, how to practically realize these basic signal processing functions inmicroscale/nanoscale has been rarely studied.The signal processing functions realized in existing MC works were performed over electricalsignals using electrical devices. In [4], [11]–[14], the transmitted bit sequence was modulatedover the concentration of signalling molecules via the on/off of an air tank [4], [11], [12],electrical spray [13], and LED controlled by Arduino microcontroller boards and laptops [14].Their high dependency over electrical signals/devices can hardly fulfil the biocompatible and non-invasive requirements of biomedical applications, such as disease diagnosis and drug delivery[15], [16]. Meanwhile, the size of electrical devices can hardly meet the requirement of intra-body healthcare applications promised by MC, where fully MC functional devices are expectedto be miniaturized into microscale/nanoscale [1].In nature, signal processing functions can be realized over molecular domain by exploiting thegene expression process, where transcription factors bind with genes to either activate or represstheir expression into proteins [17]. A gene expression process can be functioned as a buffer gate ifthe transcription factor activates the protein expression [18], and can be functioned as a NOT gateif a gene expression is repressed by the transcription factor [19]. This signal processing naturemotivates biologists to design more complex computing artificial genetic circuits to manipulatemolecular concentrations using the synthetic biology [20]. One type of artificial genetic circuitswith computing functions is the Boolean logic inspired digital logic devices. The sharp statechange between a low and high concentration is ideal for reliable state transitions and signalintegration, making digital logic particularly useful in decision-making circuits [21]. For example,the authors in [19] designed an orthogonal AND gate and coupled it to nonspecific sensors toincrease selectivity [22]. The authors in [23] constructed a simple NOR logic gate and spatiallyconfigured multiple NOR gates to produce all possible two-input gates, which have found theirutilities in biotechnological applications [24].Although the aforementioned genetic circuits has advantages in biocompatibility and minia-turization over electrical circuits, the experimental realization of genetic circuits for signalprocessing faces challenges, such as slow speed, unreliability, and non-scalability [25]. Thismotivates our initial work on chemical reactions-based microfluidic circuits [26]–[28]. Unlikegenetic circuits, chemical circuits are much easier to be controlled, and their integration withmicrofluidic devices brings advantages in low reagents consumption, rapid analysis, and highefficiency [29]. In [26], [27], we designed an MC microfluidic transceiver based on chemical reactions to successfully realize the binary concentration shift keying (BCSK) modulation and de-modulation functions. The signal processing capability of chemical reactions-based microfluidiccircuits is further exploited in [28], where we provided the designs of the AND gate, NAND gate,NOR gate, OR gate, and XOR gate, which are validated through COMSOL simulations. Echoingthe discussion in [28], one challenge in realizing signal processing functions via chemicalreactions-based microfluidic circuits design is the theoretical characterization of the logic gate,which facilitates the design parameters selection for the expected gate outputs. Although wemathematically modelled the dynamics of molecular species in microfluidic channels in [26],[27], this analysis is not scalable with the increase in the number of microfluidic circuits.Motivated by above, the main contributions of this paper are as follows: • We first present a chemical reactions-based microfluidic AND gate design, based on which,we design the microfluidic transmitter and receiver with quadruple concentration shift keying(QCSK) modulation and demodulation functionalities, to show how logic computations canprocess molecular concentrations. Compared with the BCSK transceiver in [27], the QCSKtransceiver can achieve two bits signal transmission with improved data rate. • We develop a novel mathematical framework to characterize our proposed microfluidic cir-cuits, which can be applied to analyse other new and more complicated microfluidic circuits.To do so, we first analyse the concentration and velocity changes under fluid convergenceand separation, and we derive the transfer function of a straight convection-diffusion channel.Based on these, we derive the spatial-temporal concentration distribution of a convection-diffusion-reaction channel with either a thresholding reaction or an amplifying reaction. • To evaluate our proposed microfluidic designs, we identify four elementary microfluidicblocks of the basic AND gate, and define five corresponding operators to represent theoutput concentration distribution of each elementary block. Relying on this, we not onlyderive the output concentration distribution of the proposed AND gate, but also those for ourdesigned QCSK transmitter and receiver. The functionalities of our proposed microfluidicdesigns and the corresponding theoretical results are validated via simulations performed inCOMSOL Multiphysics finite element solver.The remainder of this paper is organized as follows. In Sec. II, we provide the basic mi-crofluidic channel analysis. In Sec. III, we establish a mathematical framework to theoreticallycharacterize our proposed AND gate. In Sec. IV, we propose the designs and analysis of the
QCSK transmitter and receiver. Numerical results in Sec. V validate the proposed microfluidicdesigns and their theoretical analysis. Finally. Sec. VI concludes the paper.II. B
ASIC M ICROFLUIDIC C HANNEL A NALYSIS
With the ultimate goal of designing and analysing a microfluidic system with modulationand demodulation functionalities in this paper, the basic characteristics of fluids in microfluidicchannels must be first understood. To do so, we first analyse the concentration and velocitychanges due to the hydrodynamic fluid convergence in Sec. II-A and the fluid separation in Sec.II-B. We then present the molecular concentration distribution of a convection-diffusion channelin Sec. II-C.For a Poiseuille flow travelling along the x direction of a microfluidic channel with rectangularcross-section, the velocity profile can be obtained by solving the Navier-Stokes equation, whichis [30] v ( y, z ) = 4 h ∆ pπ ηL ∞ (cid:88) n, odd n [1 − cosh ( nπyh ) cosh ( nπw h ) ] sin ( nπzh ) , (1)where ∆ p/L denotes the pressure difference between two ends of a channel with length L , η isthe fluid dynamic viscosity, and w and h are the width and the height of the cross-section.The fluid velocity follows a “parabolic” distribution, where the closer to the channel centre,the larger the fluid velocity, resulting in the maximum velocity v max occurring at the centre ofthe channel. Alternatively, it is common to use the average velocity v eff to describe the fluidvelocity, which can be presented as v eff = Qwh . (2)In (2), Q is the volumetric flow rate, which represents the fluid volume that passes per unit time,and can be calculated as the integral of velocity contributions from each lamina using [30] Q = (cid:90) S v ( y, z )d y d z = 8 h w ∆ pπ ηL ∞ (cid:88) n =1 , , , ··· [ 1 n − hπwn tanh ( nπw h )] . (3) A. Fluid Convergence at Combining Connections
In a microfluidic circuit, fluids flowing in different channels can converge to a single flow ata combining connection, and we name this behaviour as fluid convergence for simplicity.
Inlet 2 Inlet 3 Inlet 1
Inlet n
CH1 CH2 ...
CH(n-1)
Fig. 1. A microfluidic device for fluid convergence analysis.
1) Concentration Change:
Let us consider a microfluidic device with n inlets and n − combining channels as shown in Fig. 1. We assume that species S i (1 ≤ i ≤ n ) is constantlyinjected into Inlet i with concentration C S i , average velocity v eff i , and volumetric flow rate Q i . According to the well-known analogy between hydraulic circuits and electrical circuits, thepressure drop, the flow rate, and the flow resistance in hydraulic circuits are analogous to thevoltage drop, the electrical current, and the electrical resistance in electrical circuits, respectively.Based on the Kirchhoffs Current Law , the volumetric flow rate in the i th combining channel Q CH i is the sum of coming flow rates, such that Q CH = Q + Q ,Q CH = Q CH + Q , · · · · · · ,Q CH ( n − = Q CH ( n − + Q n . (4)Therefore, the mixed concentrations of species S and S in the first combining channel are [31] C CH S = Q Q + Q C S ,C CH S = Q Q + Q C S . (5)Similarly, the mixed concentrations of species S , S , · · · , S i in the ( i − th combining channelbecome C CH ( i − S = Q + Q + ··· + Q i − ( Q + Q + ··· + Q i − )+ Q i C CH ( i − S ,C CH ( i − S = Q + Q + ··· + Q i − ( Q + Q + ··· + Q i − )+ Q i C CH ( i − S , · · · · · · ,C CH ( i − S i = Q i ( Q + Q + ··· + Q i − )+ Q i C S i . (6) Lemma 1.
For the fluid convergence from n inlets to one combining channel, the mixed con-centration of species S i can be derived as C CH( i − S i = Q i (cid:80) ij =1 Q j C S i , (7) where Q i and C S i are the volumetric flow rate and the species concentration injected into Inlet i . If all the species are injected with volumetric flow rate Q , i.e., Q = · · · = Q n = Q , species S i will be diluted to /i of its injected concentration in the ( i − th combining channel, that is C CH( i − S i = 1 i C S i . (8) Proof.
The last line of (6) can be easily reduced to (7).
Remark 1.
From (8) , we can conclude that the more inlet channels with the same volumetricflow rate injected, the lower the concentration is.2) Velocity Change:
Injecting fluids into a combining channel influence not only the speciesconcentration but also the velocity profile.
Lemma 2.
For the fluid convergence from n inlets to one combining channel, the flow rate inthe combining channel can be expressed in terms of average velocity and channel geometry as w CH( i − h CH( i − v CH( i − = w CH( i − h CH( i − v CH( i − + w i h i v eff i , (9) where v eff i , w i , and h i are the average velocity, the width, and the height of Inlet i , and v CH i eff , w CH i , and h CH i are the average velocity, the width, and the height of the i th combining channel,respectively. If all inlets and combining channels share the same geometries and the same averagevelocity v eff , the average velocity in the ( i − th combining channel becomes v CH( i − = iv eff . (10) Proof.
Based on the
Kirchhoffs Current Law and (2), we can obtain (9).
Remark 2.
It is revealed in (10) that the more inlet channels with the same volumetric flowrate injected, the larger the average velocity is.B. Fluid Separation at Bifurcation Connections
In a microfluidic circuit, a single flow can be separated into different flow streams at abifurcation connection, and we name this behaviour as fluid separation for simplicity. Let usconsider a microfluidic device with one inlet and n outlets as shown in Fig. 2(a), where a singleflow is separated into n streams travelling over n daughter channels. Assuming that species S isinjected with concentration C S and average velocity v eff , the concentration at each outlet is thesame as C S , because species S is not diluted by other species. However, the average velocityin each outlet varies for different geometries of its daughter channel. To derive the outlet averagevelocities, we establish the hydraulic circuit model in Fig. 2(b). Analogous to current division Inlet Node A Outlet nDaughter Channel n Outlet 1Daughter Channel 1 Outlet 2Daughter Channel 2 Outlet n-1Daughter Channel n-1 ... (a) A microfluidic device
Node A ... (b) Hydraulic circuit analogFig. 2. A microfluidic device for fluid separation analysis. in electric circuits, the relationship between the volumetric flow rate Q i (1 ≤ i ≤ n ) and thesupplied volumetric flow rate Q can be described by [31, Eq. (18)] Q i = R eq R i Q, (11)where R i is the hydraulic resistance of the i th daughter channel and R eq is the equivalentresistance of all daughter channels. Let us denote L D i as the length from the crosspoint NodeA in Fig. 2(b) to outlet i , and w i and h i as the geometry width and height of the i th daughterchannel, R i [32, Eq. (6)] and R eq can be calculated as R i = 12 ηL D i w i h i (cid:18) − ∞ (cid:80) i =1 , , , ··· h i wπ i tan h i ( iπh i w i ) (cid:19) , (12)and R eq = 1 / (1 /R + 1 /R + · · · + 1 /R n ) . (13) Lemma 3.
For the fluid separation from one inlet to n outlets, the average velocity v eff i in the i th outlet can be derived as v eff i = 1 L D i (cid:80) ni =1 1 L Di whw i h i v eff . (14) If all daughter channels share the same geometries, (14) can be reduced to v eff i = 1 n v eff . (15) Proof.
Combining (2), (11)-(13), we can obtain (14).
Remark 3.
It is indicated from (15) that fluid separation results in a reduction of averagevelocity by n times. C. Convection-Diffusion Channel
When a flow containing species S i enters the dispersion regime [33], we can describe thespatial-temporal concentration distribution of species S i using a 1D convection-diffusion equationas ∂C S i ( x, t ) ∂t = D eff ∂ C S i ( x, t ) ∂x − v eff ∂C S i ( x, t ) ∂x , (16)where D eff is the Taylor-Aris effective diffusion coefficient. For a microfluidic channel withrectangular-shaped cross-section whose height is h and width is w , D eff can be calculated as[34] D eff = 1 + 8 . v eff h w D ( h + 2 . hw + w ) , (17)where D is the molecular diffusion coefficient.Although the solution of (16) with a rectangular input has been derived in [27], its complexexpression does not allow the cascaded channels to be mathematically solvable in closed-form.This motivates us to derive the transfer function of a microfluidic channel so that the output ofa microfluidic circuit can be written as the convolution of an input and a cascade of the transferfunction of each channel. We solve the transfer function in the following theorem. Theorem 1.
The transfer function of a straight convection-diffusion channel is derived as H ( x, t ) = 12 π (cid:90) ∞ [ e − jωt (cid:102) C S i ( x, ω ) + e jωt (cid:102) C S i ( x, ω )]d w, (18) where (cid:102) C S i ( x, s ) = exp (cid:2) v eff x/ D eff − (cid:113) x ( v + 4 jwD eff ) / D eff 2 (cid:3) . Proof.
See the Appendix A.From
Theorem 1 , the outlet concentration of species S i can be expressed as C S i ( x, t ) = C S i ( t ) ∗ H ( x, t ) , (19)where C S i ( t ) is the input concentration of species S i at channel inlet and “ ∗ ” denotes theconvolution operator. III. AND L OGIC G ATE D ESIGN AND A NALYSIS
In this section, we present the design of the AND logic gate to demonstrate the logic com-putation ability of microfluidic circuits. The chemical reactions used in the AND gate can becategorized into two forms: 1) the thresholding reaction S i + S j → S k , and 2) the amplifyingreaction S i + Amp → S i + O . To characterize the outlet concentration of our designed AND gate,we first study the concentration distribution of the reaction channel with either a thresholding Fig. 3. The chemical reactions-based microfluidic AND logic gate. To distinguish convection-diffusion channels and convection-diffusion-reaction channels, the latter are filled with grey-gradient colour. reaction, or an amplifying reaction. Relying on the analysis in Sec. II, we then define and modelfour elementary blocks in order to theoretically characterize the AND gate.The proposed
AND gate design is presented in Fig. 3. As shown in Fig. 3, the proposed ANDgate consists of the input species I and I , and the output species O . Throughout this paper, weuse non-zero concentration to represent HIGH state (bit-1), and zero concentration to representLOW state (bit-0). The two input species I and I are first converted to an intermediate species N , so the concentration of species N will be HIGH state if either I or I is HIGH. Then,species N is further depleted by species T hL to extract the interval where both the inputs areHIGH. Finally, the remaining species N catalyses the conversion of species Amp to outputspecies O . A. Channel Model with a Thresholding Reaction
For a straight microfluidic channel with the thresholding reaction S i + S j → S k , the spatial-temporal concentration distributions of reactant and product can be expressed by convection-diffusion-reaction equations, which are expressed as ∂C S i ( x, t ) ∂t = D eff ∂ C S i ( x, t ) ∂x − v eff ∂C S i ( x, t ) ∂x − kC S i ( x, t ) C S j ( x, t ) , (20) ∂C S k ( x, t ) ∂t = D eff ∂ C S k ( x, t ) ∂x − v eff ∂C S k ( x, t ) ∂x + kC S i ( x, t ) C S j ( x, t ) , (21)where k is the rate constant. Compared with the convection-diffusion equation in (16), thenewly introduced reaction term is fully coupled with convection and diffusion process, whichcomplicates the resolution of (20) and (21). A strategy to tackle this coupling is to apply the“operator splitting” method. It first separates an original differential equation into several parts,then separately computes the solution of each part, and finally combines these separate solutionsto form a solution for the original equation. As we already derived the transfer function of a straight convection-diffusion channel in Theorem 1 , this motivates us to separate a convection-diffusion-reaction equation into a reaction term and a convection-diffusion term. This separationcan be achieved via 1) assuming the reactants are added into a virtual reactor, and the unconsumedreactants and generated product flow into a convection-diffusion channel as soon as the reactionstops; and 2) treating the solution of the reaction term as the initial input for the convection-diffusion part.With species S i and S j continuously flowing into a channel, we regard that S i and S j arecontinuously added into a virtual reactor, where the continuous reactant supply is a superpositionof reactant addition with constants at different times. To solve this, we consider the followingtwo scenarios: • Scenario 1 : species S i and S j are only added at t = 0 with concentration C S i and C S j , • Scenario 2 : species S i and S j are added continuously with concentration C S i ( t ) and C S j ( t ) .We first derive concentration changes of reactants and product for Scenario 1 , which will thenbe applied in
Scenario 2 to derive the solutions of the separated reaction term. Scenario 1 : Let c ( t ) denote the consumed concentration of reactant S i /S j during thereaction. It is noted that c ( t ) can also represent the concentration of generated species S k duringthe reaction due to a one-to-one stoichiometric relation between reactants and product. Theremaining concentrations of species S i and S j can be expressed as C S i ( t ) = C S i − c ( t ) , (22a) C S j ( t ) = C S j − c ( t ) . (22b)Then, the reaction equation can be expressed as [35, Eq. (9.13)] d[ C S i − c ( t )]d t = − k [ C S i − c ( t )][ C S j − c ( t )] . (23)After rearrangement, (23) becomes d c ( t )[ C S i − c ( t )][ C S j − c ( t )] = k d t. (24)By taking the integral of the two sides of (24), we yield c ( t ) = C Si C Sj exp [( C Sj − C Si ) kt ] − C Si C Sj C Sj exp [( C Sj − C Si ) kt ] − C Si , C S i ≤ C S j , C Si C Sj exp [( C Si − C Sj ) kt ] − C Si C Sj C Si exp [( C Si − C Sj ) kt ] − C Sj , C S i ≥ C S j . (25) Algorithm 1:
The Calculation of Remaining Concentrations of Species S i and S j Input:
The input concentrations C S i ( t ) and C S j ( t ) . The calculation time interval [0 , T ] .The time step ∆ t . Initialization of C S i = C S i ,a and C S j = C S j ,a . for n ← , (cid:98) T / ∆ t (cid:99) do Calculate the consumed concentration c n − during [( n − t, n ∆ t ] according to (25)by interchanging C S i → C n − S i and C S j → C n − S j . Update the remaining concentration C nS i ,r = C n − S i − c n − and C nS j ,r = C n − S j − c n − . Update the initial concentration C nS i = C nS i ,r + C nS i ,a and C nS j = C nS j ,r + C nS j ,a for [ n ∆ t, ( n + 1)∆ t ] . endRemark 4. It can be observed from (25) that c ( t ) is proportional to the rate constant k . Thehigher the rate constant, the faster a reactant is consumed and decreased to zero. Lemma 4.
When the reaction rate k → ∞ , the consumed concentration c ( t ) of reactant S i /S j with thresholding reaction S i + S j → S k in the microfluidic channel can be derived as lim k →∞ c ( t ) = ϕ ( C S i , C S j ) , (26) where C S i and C S j are the injected concentrations of specie S i and S j , and ϕ ( · , · ) is definedas ϕ ( x, y ) = min { x, y } . (27) Proof.
With k → ∞ , (25) can be easily reduced to (26). Scenario 2 : We consider the continuous injection of species S i and S j with concentration C S i ( t ) and C S j ( t ) . Scenario 2 can be regarded as a superposition of
Scenario 1 in timedomain. To apply the analysis of
Scenario 1 , we first discretize the reaction process into manytime intervals with the step ∆ t . Thus, the added concentration of species S i can be denoted as C nS i ,a = C S i ( n ∆ t ) ( n ≥ , where the subscript a refers to addition. We also denote C nS i and C nS i ,r as the initial and the remaining concentrations of S i at t = n ∆ t , respectively. The samenotations are also applied to species S j .We propose Algorithm 1 to numerically calculate the remaining concentrations of S i and S j after reaction S i + S j → S k . Algorithm 1 describes that for any time interval [ n ∆ t, ( n +1)∆ t ] , the consumed concentration can be calculated according to (23), but with different initial Fig. 4. Illustration of reaction S i + S j → S k in a bottle-shaped virtual reactor and a microfluidic channel. The two reactantsare marked with different colours. concentrations C nS i . This is due to the fact that the initial concentration at any time intervalis not only influenced by the newly added concentration, but also the incompletely consumedconcentration that added in previous intervals. For instance, the initial concentration C S i forthe time interval [∆ t, t ] is the sum of the newly added concentration C S i ,a and the remainingunreacted concentration C S i ,r that added at t = 0 .It is noted that the value of the rate constant k influences the accuracy of the approximation.The smaller the k , the larger volume of reactants remain. The unconsumed reactants accumulatein reactor and would participate into the reaction of the next time interval, which introducescorrelation between different time intervals. By contrast, this correlation does not exist in practicalscenario. As shown in Fig. 4, for time interval [∆ t, t ] , the flowing fluid carries remainingreactants added at t = 0 and t = ∆ t forward, which prevents them from interacting with eachother. Therefore, in the virtual reactor, we should make the rate constant approach infinity toensure reaction is always completed inside any time interval, thus eliminating the correlation.With ∆ t → , the remaining concentrations of S i and S j calculated in Algorithm 1 reduce to C S i ,r ( t ) = C S i ( t ) − ϕ [ C S i ( t ) , C S j ( t )] , (28a) C S j ,r ( t ) = C S j ( t ) − ϕ [ C S i ( t ) , C S j ( t )] , (28b)where ϕ ( · , · ) is given in (27).We derive the output concentration distributions of species S i , S j and S k in the followinglemma. Lemma 5.
For a straight reaction channel with the thresholding reaction S i + S j → S k , the output concentration distributions of species S i , S j , and S k can be derived as C S i ( x, t ) = C S i ,r ( t ) ∗ H ( x, t ) , (29a) C S j ( x, t ) = C S j ,r ( t ) ∗ H ( x, t ) , (29b) C S k ( x, t ) = ϕ [ C S i ( t ) , C S j ( t )] ∗ H ( x, t ) , (29c) where C S i ,r ( t ) , C S j ,r ( t ) , H ( x, t ) , and ϕ ( · , · ) are given in (28a) , (28b) , (18) , and (27) , respectively.Proof. Recall that we separate a convection-diffusion-reaction equation into a reaction part anda convection-diffusion part, we consider the remaining concentrations of S i in (28a) and S j in(28b) as inputs to a straight convection-diffusion channel. According to (19), we can obtain (29a)and (29b). The derivation of (29c) can see Appendix B. B. Channel Model with an Amplifying Reaction
Lemma 6.
For a straight reaction channel with the amplifying reaction S i + Amp → S i + O ,the output concentration distribution of species O can be derived as C O ( x, t ) = (cid:2) C Amp ( t ) · { C Si ( t ) > } (cid:3) ∗ H ( x, t ) , (30) where C Amp ( t ) and C S i ( t ) are the injected concentrations of species Amp and S i , u ( t ) is theHeaviside step function, {·} is the indicator function that represents the value 1 if the statementis true, and zero otherwise.Proof. To analyse a straight microfluidic channel with amplifying reaction S i + Amp → S i + O , we also separate it into a reaction term and a convection-diffusion term. For the reactionterm, as species O is only produced in the presence of S i and the concentration equals theinjected concentration of species Amp [36], the reaction solution can be expressed as C Amp ( t ) · { C Si ( t ) > } . Taking the reaction solution as the initial input for a convection-diffusion channel,the concentration of product O can be derived as (30). C. Elementary Blocks
Relying on the analysis of fluid convergence in
Lemma 1 and 2 , convection-diffusion channelin
Theorem 1 , and convection-diffusion-reaction channel in
Lemma 5 and 6 , we focus on theanalysis of four elementary blocks in our designed AND gate (Fig. 3) in Table I. Meanwhile,we define five typical operators for the four elementary blocks, aiming at simplifying the outputexpression of a microfluidic circuit. As shown in Table I, the operator T [ · ] represents the output TABLE IF
OUR ELEMENTARY BLOCKS . Operator Elementary Block Operator Output T [ C S i ( t ) , n ] Eq. (31) C S i ( t ) : The output of aconvection-diffusion channelwith length L T . G [ C S i ( t ) , C S j ( t ) , n ] Eq. (33) C S k ( t ) : The concentration ofproduct S k with S i + S j → S k . R [ C S i ( t ) , C S j ( t ) , n ] Eq. (34) C S i ( t ) : The remaining concentrationof S i with S i + S j → S k . A [ C S i ( t ) , C Amp ( t ) , n ] Eq. (35) C O ( t ) : The concentration ofproduct O with S i + Amp → S i + O . F [ C S i ( t ) , C S j ( t ) , C Amp ( t ) , n ] Eq. (36) C O ( t ) : The concentration ofproduct O with S i + S j → S k and S i + Amp → S i + O . of a convection-diffusion channel with length L T , and can be expressed as T [ C S i ( t ) , n ] (cid:44) C S i ( t ) ∗ H n ( L T , t ) , (31)where the subscript n of H n indicates that the average velocity in the channel is nv eff .For the block with thresholding reaction S i + S j → S k , species S i with concentration C S i ( t ) and velocity ( n − v eff , and species S j with concentration C S j ( t ) and velocity v eff , are injectedto the channel with length L C . The convergence of two subchannels with species S i and S j willresult in a concentration dilution, and the diluted concentrations of S i and S j are ( n − C S i ( t ) /n and C S j ( t ) /n following (7) in Lemma 1 , respectively. Meanwhile, the average velocity willincrease to nv eff following (9) in Lemma 2 . Then, species will flow to a buffer channel beforethe convection-diffusion-reaction channel filled with grey-gradient colour. The buffer channelallows the reactants to be well mixed before a reaction, and the reactant mixing along the radialdirection only relies on diffusion. To achieve a fully diffusional mixing, the minimum bufferlength L B can be estimated as L B = w + h D v eff , (32)where w + h D is the maximum time to travel over the radial direction. We define operator G [ · ] todescribe the concentration of product S k , and according to (29c), G [ · ] can be expressed as G [ C S i ( t ) , C S j ( t ) , n ] (cid:44) ϕ [( n − C S i ( t ) /n, C S j ( t ) /n ] ∗ H n ( nL B + L C , t ) ∗ H n ( L R , t ) . (33) For the same reaction, we define operator R [ · ] to characterize the residual concentration of S i .According to (29a), R [ · ] can be expressed as R [ C S i ( t ) , C S j ( t ) , n ] (cid:44) (cid:2) ( n − C S i ( t ) /n − ϕ [( n − C S i ( t ) /n, C S j ( t ) /n ] (cid:3) ∗ H n ( nL B + L C , t ) ∗ H n ( L R , t ) . (34)For the amplifying reaction S i + Amp → S i + O , operator A [ · ] describes the concentration ofproduct O , and can be expressed using Lemma 6 as A [ C S i ( t ) , C Amp ( t ) , n ] (cid:44) (cid:2) [ C Amp ( t ) /n ∗ H n ( nL B + L C , t )] · { [( n − C Si ( t ) /n ∗ H n ( nL B + L C ,t )] > } (cid:3) ∗ H n ( L R , t ) . (35)As seen in the AND gate design in Fig. 3, a threshold reaction is cascaded with an amplifyingreaction; thus, we define operator F [ · ] as a combination of operators R [ · ] and A [ · ] , whichrepresents the concentration of product O with S i + S j → S k and S i + Amp → S i + O as F [ C S i ( t ) , C S j ( t ) , C Amp ( t ) , n ] (cid:44) A (cid:2) R [ C S i ( t ) , C S i ( t ) , n ] , C Amp ( t ) , n + 1 (cid:3) . (36) D. AND Logic Gate Analysis
We denote the concentrations of input species I and I as C I ( t ) and C I ( t ) . Remindthat we use non-zero concentration to represent HIGH state (bit-1), and zero concentration torepresent LOW state (bit-0). Therefore, we assume that at any time t , C I ( t ) and C I ( t ) eitherare HIGH concentration C or LOW concentration . Species M , T hL , and
Amp are injectedcontinuously; thus, their initial concentrations follow C M ( t ) = C M u ( t ) , C T hL ( t ) = C T hL u ( t ) ,and C Amp ( t ) = C Amp u ( t ) . For simplicity, all reactants are injected using a same averagevelocity v eff . Theorem 2.
The concentration distribution of the product species O in our designed AND gatein Fig. 3 can be derived as C O ( x , t ) = F (cid:8) T [ C N ( x , t ) , , T [ C T hL ( t ) , , T [ C Amp ( t ) , , (cid:9) , (37) where C N ( x , t ) = 12 (cid:8) G (cid:2) T [ C I ( t ) , , T [ C M ( t ) , , (cid:3) + G (cid:2) T [ C I ( t ) , , T [ C M ( t ) , , (cid:3)(cid:9) ∗ H ( L A , t ) . (38) In (37) and (38) , operators T [ · ] , G [ · ] , F [ · ] are defined in Table I, L A is the travelling distanceof the laminar located at the centre channel from x to x in Fig. 3. Fig. 5. The concentration of species N before and after the reaction N + T hL → W in the design AND gate. Proof.
As shown in Fig. 3, the species N generated by two inputs merge with each other afterreactions I M → N and I M → N at position x = x . With L A , the concentrationof species N at x = x can be derived as (38). Then, species N travels over a convection-diffusion channel and enters the elementary block F [ · ] consisting of reactions N + T hL → W and N + Amp → N + O to produce the gate output O . According to the definition of F [ · ] in(36), the concentration of species O at location x = x can be derived as (37).For the thresholding reaction N + T hL → W in Fig. 3, C T hL directly determines the gatefunction. We derive the constraint for C T hL in the following lemma. Lemma 7.
To ensure that our designed AND gate exhibits desired behaviour, the concentrationof species
T hL needs to satisfy lim t →∞ ϕ ( C , C M ) T [ u ( t ) , ∗ q ( t ) T [ u ( t ) , ∗ H (5 L B + L C , t ) < C T hL < t →∞ ϕ ( C , C M ) T [ u ( t ) , ∗ q ( t ) T [ u ( t ) , ∗ H (5 L B + L C , t ) , (39) where C is the HIGH concentration of input species I and I , C M is the injected concentrationof species M , q ( t ) = H ( L A , t ) ∗ H ( L T , t ) ∗ H (5 L B + L C , t ) , H ( x, t ) is the transfer functionof a convection-diffusion channel derived in (18) , and T [ · ] is defined in (31) .Proof. Let C and C T hL denote the steady-state concentrations of species N and T hL at location x = x , respectively. Fig. 5 plots the concentration of species N before and after reaction N + T hL → W . When only one input is HIGH, the steady-state concentration C can beexpressed as C = lim t →∞ · · G (cid:8) T [ C u ( t ) , , T [ C M ( t ) , , (cid:9) ∗ q ( t ) . (40)When both inputs are HIGH, the steady-state concentration becomes C . For species T hL , itssteady-state concentration C T hL at x = x can be expressed as C T hL = lim t →∞ T [ C T hL u ( t ) , ∗ H (5 L B + L C , t ) . (41)As shown in Fig. 5, the blue region represents that both two inputs are HIGH, and the yellowregion represents that only one input is HIGH. The relationship between C and C T hL has threecases: • C T hL < C : After reaction, the remaining concentration of species N contains the regionwhere one or both the inputs are HIGH. • C < C T hL < C : After reaction, the remaining concentration of species N only containsthe region where both two inputs are HIGH. • C < C T hL : After reaction, species N is completely depleted.Therefore, to capture the region where both the inputs are HIGH, the concentration of species T hL needs to satisfy the condition C < C T hL < C . Combined with (40) and (41), we canderive (39). IV. M ICROFLUIDIC
QCSK T
RANSMITTER AND R ECEIVER
In this section, we present the microfluidic designs to show how logic computations can processmolecular concentration so as to achieve QCSK modulation and demodulation. Meanwhile,we also theoretically characterize the output concentration distributions of the proposed QCSKtransmitter and receiver.
A. QCSK Transmitter1) QCSK Transmitter Design:
QCSK modulation represents two digital inputs as four con-centration levels of an output signal, which is analogous to the Amplitude Shift Keying (ASK)modulation in wireless communication [3]. A challenge of implementing a QCSK MC transmitteris how to control the output concentration via four different input combinations, i.e. , “00”, “01”,“10”, and “11”. We solve this challenge by borrowing the idea of an electric 2:4 decoder. Inelectric field, a 2:4 decoder, which has 2 inputs and 4 outputs, selects exactly one of its outputsaccording to the input combination. Fig. 6 presents the truth table and an implementation for anelectric 2:4 decoder, where four AND gates receive the HIGH or the LOW of input species I and I .Inspired by the electric 2:4 decoder, we propose a chemical reactions-based microfluidic 2:4decoder to realize QCSK modulation as Fig. 7. The proposed microfluidic device is made upof four microfluidic units corresponding to four different concentration outputs. For ease ofreference, these four units are named as unit 4, unit 3, unit 2, and unit 1 from top to bottom. Fig. 6. The truth table and implementation of an electric 2:4 decoder.
Fig. 7. The chemical reactions-based microfluidic 2:4 decoder.
Analogous to the electric 2:4 decoder in Fig. 6, the AND gate in each unit either takes I and I or their complementary species P and P as its inputs. It is noted that species P and P aresupplied continuously with a HIGH state so that after reactions I P → W and I P → W , the remaining concentrations of species P and P , i.e. , C P /P ( x , t ) , can represent thecomplementary states of species I and I , thus achieving the NOT gate. Unlike an electric 2:4decoder that an identical voltage level is produced no matter which unit is selected, the proposedchemical 2:4 decoder will output different concentration levels. As each unit output C iO ( t ) isinfluenced by C iAmp in an amplifying reaction, the concentration variation of transmitted signalsis represented via different concentrations of injected species Amp as C iAmp ( t ) = C iAmp u ( t ) ( ≤ i ≤ ) for different units. Here, we set C Amp > C Amp > C Amp > C Amp to ensure max { C O ( t ) } > max { C O ( t ) } > max { C O ( t ) } > max { C O ( t ) } .
2) QCSK Transmitter Analysis:
The objective of the following analysis is to derive thetransmitter output C iO ( t ) of the design in Fig. 7. We first derive the inputs of AND gates, i.e. , the concentrations of I /I and P /P at location x = x . When input species I and I directly flow into an AND gate, their concentrations at location x = x can be expressed as C I /I ( x , t ) = (cid:2) T [ C I /I ( t ) , ∗ H ( L C + L B + L R , t ) (cid:3) / , (42)where C I /I ( t ) is the concentration of input species I /I , operator T [ · ] is defined in TableI. When the complementary species P and P flow into an AND gate, their concentrations atlocation x = x can be expressed as C P /P ( x , t ) = R (cid:8) T [ C P /P ( t ) , , T [ C I /I ( t ) , , (cid:9) , (43)where C P /P ( t ) is the concentration of species P /P , operator R [ · ] is also defined in TableI. With the derived AND gate inputs C I /I ( x , t ) in (42) and C P /P ( x , t ) in (43), the trans-mitter output C iO ( t ) can be expressed using Theorem 2 by interchanging the parameters • in (38) via: T [ C I ( t ) , → C I /I ( x , t ) if an AND gate input is I /I , T [ C I ( t ) , → C P /P ( x , t ) if an AND gate input is P /P , n = 2 → n = 3 , H ( L A , t ) → H ( L A , t ) ; • in (37) via: T [ C N ( x , t ) , → T [ C N ( x , t ) , , C Amp ( t ) → C iAmp ( t ) , n = 5 → n = 7 . B. QCSK Receiver1) QCSK Receiver Design:
From the communication perspective, a corresponding microflu-idic receiver is required to distinguish different concentration levels of C iO ( x , t ) from differentinput combinations to achieve QCSK demodulation. In the following, we simplify C iO ( x , t ) using C iO ( t ) to represent a selected output of the proposed QCSK transmitter, and C O ( t ) torepresent the general receiver input. In this paper, we consider the setup where the QCSKtransmitter output C iO ( t ) directly flows into the QCSK receiver; therefore, the receiver input C O ( t ) ∈ (cid:2) , max { C O ( t ) } (cid:3) . We also denote C Y ( t ) and C Y ( t ) as the final demodulated concen-tration signals, which correspond to the transmitter concentration inputs C I ( t ) and C I ( t ) ,respectively.To detect four concentration levels at the output of our proposed QCSK transmitter, we firstdesign three detection microfluidic units in Fig. 8 to serve as a front-end processing module forthe QCSK receiver, where each detection unit follows the receiver design in our initial work [27],with the capability of generating a rectangular output if the maximum concentration of a received Fig. 8. Three detection units [27] serve as a front-end processing modules. Each channel is labelled with a channel numberto denote channel length as L number . By setting max { C iO ( t ) } < C iT < max { C i +1 O ( t ) } , the front-end processing module candistinguish four concentration regions. TABLE IIT HE RELATION BETWEEN THE RECEIVER INPUT C O ( t ) , FRONT - END MODULE OUTPUT BINARY SIGNAL B , AND RECEIVEROUTPUT BINARY SIGNAL Y . max { C O ( t ) } B B B Y Y [0 , C T ] 0 0 0 0 0[ C T , C T ] 0 0 1 0 1[ C T , C T ] 0 1 1 1 0[ C T , ∞ ) 1 1 1 1 1 signal exceeds a predefined threshold. As shown in Fig. 8, the only difference among the threedetection units is the injected concentration C iT ( t ) = C iT u ( t ) (1 ≤ i ≤ of thresholdingreactant T . By setting max { C iO ( t ) } < C iT < max { C i +1 O ( t ) } , the concentration region of C O ( t ) can be identified for three-bit binary signals B B B as shown in Fig. 8. For instance, if max { C O ( t ) } > C T , all detection units will output a HIGH state with B B B = 111 .It is noted that the three detection units in Fig. 8 can only demodulate C O ( t ) to threeconcentration signals C B ( t ) , C B ( t ) , and C B ( t ) instead of C Y ( t ) and C Y ( t ) , which means extrasignal processing units are required. Consider the outputs of front-end module C B ( t ) exhibita rectangular concentration profile and its digital characteristic is ideal to perform logic com-putations [19], this motivates us to design logic circuits to transform C iB ( t ) to desirable output C Y ( t ) and C Y ( t ) . To inspire the design for this signal transformation, we present the relationshipbetween the binary signal B i ( ≤ i ≤ ) and the binary signal Y j ( j = 1 , ) in the truth tableof Table II. Based on Table II, we express the Boolean equations [37] for Y and Y as Front-end processing module Back-end processing module (a) Microfluidic channels to calculate C Y Front-end processing module Back-end processing module (b) Microfluidic channels to calculate C Y Fig. 9. The microfluidic QCSK receiver design. Each channel is labelled with a channel number to denote channel length as L number . The enlargement of (a) and (b) can see Appendix C and D, respectively. Y = ¯ B B B + B B B = B B , (44)and Y = ¯ B ¯ B B + B B B = B ( B (cid:12) B ) , (45)where ¯ B is the complementary form of B , B B represents the AND operation of B and B , and (cid:12) is the Exclusive NOR (XNOR) operation. Inspired by these Boolean relationshipsbetween species B and receiver output species Y in (44) and (45), we connect the front-endmodule with an AND gate to compute C Y ( t ) as shown in Fig 9(a), as well as a NXOR gateand an AND gate to calculate C Y ( t ) as shown in Fig. 9(b).
2) QCSK Receiver Analysis:
To theoretically characterize receiver outputs C Y ( t ) and C Y ( t ) ,we denote C [ · ] ( t ) as the concentration of any injected species [ · ] , and L i as the length of themicrofluidic channel with number i . Moreover, we assume that all types of molecules are injectedwith average velocity v eff . In the following, we first derive the front-end processing output C iB ( t ) in Fig. 8, and then derive the QCSK receiver outputs C Y ( t ) and C Y ( t ) in Fig. 9. In addition, the location and channel number are in bold in the following so that readers can easily followour derivation. C iB ( t ) Derivation:
As shown in Fig. 8, each detection unit in the front-end processing moduleconsists of a thresholding reaction O + T → W and an amplifying reaction O + A → O + B ,and the output can be expressed using the operator F [ · ] defined in Table I as C iB ( x , t ) = F (cid:8) T [ C O ( t ) , , T [ C iT ( t ) , , T [ C A ( t ) , , (cid:9) , (46)where C O ( t ) is the receiver input concentration. C Y ( t ) Derivation:
As shown in in Fig. 9(a), C B ( t ) and C B ( t ) flow into an AND gate toproduce C Y ( t ) . At x = x , C Y ( t ) can be derived as C Y ( x , t ) = F (cid:8) T [ 12 (cid:88) j =1 C jB ( x , t ) ∗ H ( 2 L + L + h , t ) , T [ C T ( t ) , , T [ C A ( t ) , , (cid:9) , (47)where / represents the dilution of C / B ( x , t ) by C / B ( x , t ) , H n ( x, t ) is given in Theorem 1 with n indicating that the average velocity is nv eff , and the operator T [ · ] is defined in Table I. C Y ( t ) Derivation:
As shown in Fig. 9(b), an XNOR gate and an AND gate are linked tothe front-end processing module to produce C Y ( t ) . • XNOR Gate Analysis : Relying on the fluid separation analysis in
Lemma 3 , at x = x , C jB ( x , t ) ( j = 2 , is equally separated from channel to channels due to thesymmetrical microfluidic design from x to x in Fig. 9(b), resulting in a velocity reductionfrom v eff in channel with O + A → O + B to . v eff in channels . In channels , theconfluence of C B ( x , t ) and C B ( x , t ) occurs, and then is diluted by species T injectedat x . Subsequently, the outer fluid performs reaction B + T → W to capture the regionwhere both C B ( x , t ) and C B ( x , t ) are HIGH as the second case in Fig. 5, while the innerfluid flows forward without this reaction. At x = x , the amplifying products R and R after reactions B + A → B + R and B + A → B + R can be expressed as C R ( x , t ) = A (cid:8) T [ 12 (cid:88) j =2 C jB ( x , t ) ∗ H . ( 3 L + L + 2 L + L + h + 2 w , t ) , ∗ H ( L + L + L , t ) (cid:124) (cid:123)(cid:122) (cid:125) C Inner B ( x ,t ) , T [ C A ( t ) , , (cid:9) , (48) and C R ( x , t ) = F (cid:8) T [ 12 (cid:88) j =2 C jB ( x , t ) ∗ H . ( 3 L + L + 2 L + L + h + 2 w , t ) , (cid:124) (cid:123)(cid:122) (cid:125) C Outer B ( x ,t ) , T [ C T ( t ) , , T [ C A ( t ) , , (cid:9) , (49)where the superscript “Inner” and “Outer” represent the outer and inter fluids from x to x , and / in (48) represents the dilution of species B by species T .After reaction R R → W , the remaining species R at x = x will be HIGH wheneither C B ( x , t ) or C B ( x , t ) is HIGH, thus achieving an XOR gate. Relying on (29a) in Lemma 5 , the remaining concentration of species R is derived as C R ( x , t ) = 12 (cid:8) C R ( x , t ) − ϕ [ C R ( x , t ) , C R ( x , t )] (cid:9) ∗ H ( L + L + 2 w , t ) ∗ H ( L , t ) ∗ H ( L , t ) , (50)where ϕ [ · , · ] is given in (27). The cascaded reaction R N OT → W functions as a NOTgate similar to the reaction I P → W in the QCSK transmitter in Fig. 7 in order toachieve the XNOR gate. At x = x , the concentration of N OT can be expressed usingthe operator R [ · ] defined in Table I as C NOT ( x , t ) = R (cid:8) T [ C NOT ( t ) , , C R ( x , t ) , (cid:9) ∗ H ( 2 L + L + h , t ) , (51)where the superscript represents the species N OT generated by C B ( x , t ) and C B ( x , t ) . • AND Gate Analysis : The calculation of receiver output C Y ( t ) also needs the participationof C B ( x , t ) . To perform the AND gate, the product species B (indicated by the red arrow)should be converted to molecular type N OT via B + V → N OT . At x = x , theconcentration of species N OT generated by C B ( x , t ) can be expressed using operator G [ · ] defined in Table I as C NOT ( x , t ) = G (cid:8) C B ( x , t ) ∗ H ( L , t ) , T [ C V ( t ) , , (cid:9) ∗ H ( 2 L + 2 L + L + h , t ) . (52)We highlight that C NOT ( x , t ) and C NOT ( x , t ) must be well synchronized. This meansthat C NOT ( x , t ) and C NOT ( x , t ) should arrive at x simultaneously, which can beachieved by ensuring the inputs C O ( t ) of three detection units have the same travellingtime from the front-end module to position x in Fig. 9(b). Finally, we can derive theQCSK receiver output C Y ( x , t ) as C Y ( x , t ) = F (cid:8) T [ 415 C NOT ( x , t ) + 1115 C NOT ( x , t ) , , T [ C T ( t ) , , T [ C A ( t ) , , (cid:9) , (53)where / represents the dilution of C NOT ( x , t ) by C NOT ( x , t ) , while / representsthe dilution of C NOT ( x , t ) by C NOT ( x , t ) .V. P ERFORMANCE E VALUATION
In this section, we implement our proposed microfluidic AND gate, QCSK transmitter, andQCSK receiver design in Fig. 3, Fig. 7, and Fig. 9 using COMSOL Multiphysics, which arethen used to validate our corresponding theoretical analysis. The transfer function H ( x, t ) givenin Theorem 1 is computed in Matlab using quadgk . As quadgk is only an approximation of H ( x, t ) , the computed results may fluctuate around their steady values. If a computed value isslightly larger than steady value , this can induce an instant change on the output value of theindicator function in (30) from to , which would further lead to a generation of output signalsin undesired regions after an amplifying reaction. To avoid this phenomenon, we modify thestatement of an indicator function C S i ( t ) > as C S i ( t ) > max { C S i ( t ) } . By doing so, thewidth of a rectangular output is expected to be smaller than that of the corresponding simulationresult. In COMSOL simulations, unless otherwise stated, we set v eff = 0 . cm/s, D eff = 10 − m /s, w = 20 µ m, h = 10 µ m, k = 400 m /(mol · s). Furthermore, we use Ana. and Sim. to abbreviateAnalytical and Simulation in all figures. A. AND Logic Gate
Fig. 10 presents the COMSOL simulation results of the AND logic gate design depicted inFig. 3. We set the parameters: C I ( t ) = 8[ u ( t − − u ( t − , C I ( t ) = 8[ u ( t − − u ( t − , C M ( t ) = 8 u ( t ) , C Amp ( t ) = 12 u ( t ) , L T = 80 µ m, L C = 20 µ m, L R = 500 µ m, L A = 120 µ m.For the injected concentration of species T hL , we consider three cases: C T hL ( t ) = 5 u ( t ) , C T hL ( t ) = 10 u ( t ) , C T hL ( t ) = 20 u ( t ) , with the aim to examine its impact on the gate behaviour.Fig. 10(a) plots the concentrations of species T hL before reaction N + T hL → W in Fig.3. We observe that the simulated concentration points agree with the analytical concentrationcurves, thus demonstrating the correctness of our analysis of convection-diffusion in Theorem 1 and convection-diffusion-reaction channels in
Lemma 5 . For the three different injected concen-trations, species
T hL is nearly diluted to one-fifth of its injected concentration due to that species
T hL enters the microfluidic device via the fifth inlet, which validates the concentration analysisfor fluid convergence in
Lemma 1 . Moreover, we also plot the concentration constraint C in (a) The concentration of species N and T hL at x = x . (b) The normalized concentrations of input species I , I , and output species O .Fig. 10. The evaluation of an AND logic gate. (40) for species T hL using black dash lines. For the curves with C T hL ( t ) = 5 u ( t ) or u ( t ) , C T hL does not satisfy the concentration constraint in Lemma 7 ; as expected, the microfluidicdevice fails to achieve the AND function, which is demonstrated in Fig. 10(b). Fig. 10(b) plotsthe normalized inputs and the final output product O in (37). Only for C T hL ( t ) = 10 u ( t ) , thewidth of species O equals the width where both input species I and I are HIGH, demonstratingthe desirable behaviour of an AND gate. Furthermore, due to the modification of the indicatorfunction set, we can see the width of (37) is smaller than that of the simulation results. B. QCSK Transmitter
Fig. 11 plots the outputs of the proposed microfluidic QCSK transmitter design in Fig. 7and their analytical values C iO ( t ) in Sec. IV-A2. Species I and I are injected with either u ( t − − u ( t − representing bit or u ( t ) representing bit . For other molecular types, theirinjected concentrations are set as: C P ( t ) = C P ( t ) = 12[ u ( t − − u ( t − , C M ( t ) = 12 u ( t ) , C T hL ( t ) = 16 u ( t ) , C Amp ( t ) = 24 u ( t ) , C Amp ( t ) = 16 u ( t ) , C Amp ( t ) = 8 u ( t ) , and C Amp ( t ) =0 . The buffer channels are configured with L B = 100 µ m, L B = 150 µ m, L B = 350 µ m, and L B = 400 µ m.As shown in Fig. 11, for any input combination, only one unit outputs a HIGH signal exceptfrom the case where both I and I are LOW due to C Amp ( t ) = 0 . Moreover, the analyticalcurves always capture the simulation points, which again demonstrates the effectiveness of ourtheoretical analysis C iO ( t ) in Sec. IV-A2. As species Amp is supplied with different injectedconcentrations for each unit, we see that the selected unit reaches different concentration levels, (a) Unit 4 (b) Unit 3(c) Unit 2 (d) Unit 1Fig. 11. The output concentrations of the proposed microfluidic QCSK transmitter.TABLE IIIT HE PARAMETERS OF THE
QCSK
RECEIVER . Molecular Type Concentration (mol/m ) Molecular Type Concentration (mol/m ) A u ( t ) T u ( t ) A u ( t ) T u ( t ) A u ( t ) T u ( t ) A u ( t ) NOT u ( t ) A u ( t ) V u ( t ) TABLE IVT
HE GEOMETRY OF THE
QCSK
RECEIVER . Channel Number Length ( µ m) Channel Number Length ( µ m) Channel Number Length ( µ m) Channel Number Length ( µ m) proving that the proposed microfluidic QCSK transmitter successfully modulates input bits tothe concentration level of output species O . C. QCSK Receiver
To evaluate the proposed QCSK receiver design in Fig. 9, we consider four different rectangularconcentration profiles as the receiver input C O ( t ) , which is C O ( t ) = i [ u ( t − − u ( t − ≤ i ≤ . Accordingly, to distinguish these four concentration levels, the concentration ofspecies T for three units in Fig. 8 are set as: C T ( t ) = 2 . u ( t ) , C T ( t ) = 1 . u ( t ) , and C T ( t ) = 0 . u ( t ) . Other parameters and the geometry are summarized in Table III and IV.Fig. 12 plots the outputs of the proposed QCSK receiver design in Fig. 9 and the correspondinganalytical results of C Y ( t ) in (47) and C Y ( t ) in (53). First, we can see that although simulation (a) C O ( t ) = 3[ u ( t − − u ( t − (b) C O ( t ) = 2[ u ( t − − u ( t − (c) C O ( t ) = u ( t − − u ( t − (d) C O ( t ) = 0 Fig. 12. The output concentrations of the proposed microfluidic QCSK receiver. curves are not in precise agreement with analytical curves, the close match can still confirm thecorrectness of the mathematical characterization of C Y ( t ) in (47) and C Y ( t ) in (53). Second,we observe the width difference between analytical and simulation curves for C Y is largerthan that for C Y . This is because the modification of the statement of an indicator functionresults in the width difference in each amplifying reaction, the more amplifying reactions areutilized to compute C Y in Fig. 9(b), the wider the width difference is. Third, we see that theproposed receiver design can well demodulate the received signal C O ( t ) to two outputs C Y and C Y . Recall that we use non-zero concentration to represent HIGH state (bit-1), and zeroconcentration to represent LOW state (bit-0). We also observe that the relationship between themaximum concentration of the receiver input max { C O ( t ) } , the concentration of species T , andbinary signals Y and Y is in consistent with the truth table of Table II, which demonstrates theeffectiveness of our proposed design. VI. C ONCLUSION
In this paper, we considered the realization of quadruple concentration shift keying (QCSK)modulation and demodulation functionalities for molecular communication (MC) using chemical reactions-based microfluidic circuits. We first presented an AND gate design to demonstrate thelogic computation capabilities of microfluidic circuits, and then showed how to utilize logiccomputations to achieve QCSK modulation and demodulation functions. To theoretically char-acterize a microfluidic circuit, we established a general mathematical framework which is scalablewith the increase of circuit scale and can be used to analyse other new and more complicatedcircuits. We derived the output concentration distributions of the AND gate, QCSK transmitterand receiver designs. Simulation results obtained from COMSOL Multiphysics showed all theproposed microfluidic circuits responded appropriately to input signals, and closely matchedour derived analytical results. We believe that this paper not only provides a design principleand mathematical framework for microfluidic MC circuits, but also a foundation for harnessingsimple microfluidic logic gates to produce diverse and complex signal processing functions.A PPENDIX AP ROOF OF T HEOREM H ( x, t ) , we formulate the following initial and boundaryconditions for (16) C S i (0 , t ) = δ ( t ) , (54a) C S i ( x,
0) = 0 , x ≥ , (54b)and ∂C S i ( x, t ) ∂x | x = ∞ = 0 , t ≥ . (54c)where δ ( · ) is the Kronecker delta function. The Laplace Transform of (16) with respect to t is D eff ∂ (cid:102) C S i ( x, s ) ∂x − v eff ∂ (cid:102) C S i ( x, s ) ∂x − s (cid:102) C S i ( x, s ) = 0 . (55)The general solution for this second order differential equation can be expressed as (cid:102) C S i ( x, s ) = d e v eff + √ v eff +4 D eff s D eff x + d e v eff − √ v eff +4 D eff s D eff x , (56)where d and d are two constants. To determine d and d , we also apply Laplace Transformsto (54a) and (54c), which are (cid:102) C S i (0 , s ) = 1 , (57) ∂ (cid:102) C S i ( x, s ) ∂x | x = ∞ = 0 . (58)Constrained by these two conditions, we arrive at the particular solution for (55) as (cid:102) C S i ( x, s ) = e v eff − √ v eff +4 D eff s D eff x , (59) In order to obtain the transfer function, we need to calculate the inverse Laplace Transform of(59), i.e. , L − (cid:110) (cid:102) C S i ( x, s ) (cid:111) . However, L − (cid:110) (cid:102) C S i ( x, s ) (cid:111) is mathematically not solvable in close-form due to the complexity of (59). Here, we resort to the Gil-Pelaez theorem and consider L − (cid:110) (cid:102) C S i ( x, s ) (cid:111) as a probability density function whose characteristic function is (cid:102) C S i ( x, s ) .The cumulative distribution function (CDF) for L − (cid:110) (cid:102) C S i ( x, s ) (cid:111) can be expressed as F ( x, t ) = 12 − π (cid:90) ∞ e − jωt (cid:102) C S i ( x, ω ) − e jωt (cid:102) C S i ( x, ω )2 jω d w. (60)Take the derivative of F ( x, t ) with respect to t , we can arrive at (18).A PPENDIX BT HE DERIVATION OF THE C ONCENTRATION OF P RODUCT S PECIES S k IN (29c)To derive the concentration of product S k , we combine (20) and (21) and denote C ( x, t ) = C S i ( x, t ) + C S k ( x, t ) , which yields ∂C ( x, t ) ∂t = D eff ∂ C ( x, t ) ∂x − v eff ∂C ( x, t ) ∂x . (61)The sum concentration has the following initial and boundary conditions C (0 , t ) = C S i (0 , t ) , (62a) C ( x,
0) = 0 , x ≥ , (62b)and ∂C ( x, t ) ∂x | x = ∞ = 0 , t ≥ . (62c)As these conditions are the same as (54a)-(54c), we can write C ( x, t ) = C S i (0 , t ) ∗ H ( x, t ) . (63)Combined with (28a) and (29a), the concentration of product S k is C S k ( x, t ) = C ( x, t ) − C S i ( x, t )= ϕ [ C S i ( t ) , C S j ( t )] ∗ H ( x, t ) . (64) A PPENDIX CT HE E NLARGEMENT OF THE BACK - END PROCESSING MODULE IN F IG . 9(a) F r on t - e nd p r o ce ss i ng m odu l e B ac k - e nd p r o ce ss i ng m odu l e Fig. 13. Microfluidic channels for the QCSK output C Y . A PPENDIX DT HE E NLARGEMENT OF THE BACK - END PROCESSING MODULE IN F IG . 9(b) F r on t - e nd p r o ce ss i ng m odu l e B ac k - e nd p r o ce ss i ng m odu l e Fig. 14. Microfluidic channels for the QCSK output C Y . R EFERENCES [1] M. Kuscu, E. Dinc, B. A. Bilgin, H. Ramezani, and O. B. Akan, “Transmitter and Receiver Architectures for MolecularCommunications: A Survey on Physical Design With Modulation, Coding, and Detection Techniques,”
Proc. IEEE , vol.107, no. 7, pp. 1302–1341, July 2019.[2] V. Jamali, A. Ahmadzadeh, W. Wicke, A. Noel, and R. Schober, “Channel Modeling for Diffusive Molecular Communi-cationA Tutorial Review,”
Proc. IEEE , vol. 107, no. 7, pp. 1256–1301, July 2019.[3] M. S. Kuran, H. B. Yilmaz, T. Tugcu, and I. F. Akyildiz, “Modulation Techniques for Communication via Diffusion inNanonetworks,” in
Proc. IEEE ICC , June 2011, pp. 1–5.[4] B. Koo, C. Lee, H. B. Yilmaz, N. Farsad, A. Eckford, and C. Chae, “Molecular MIMO: From Theory to Prototype,”
IEEEJ. Sel. Areas Commun. , vol. 34, no. 3, pp. 600–614, March 2016.[5] H. B. Yilmaz, A. C. Heren, T. Tugcu, and C. Chae, “Three-Dimensional Channel Characteristics for MolecularCommunications With an Absorbing Receiver,”
IEEE Commun. Lett. , vol. 18, no. 6, pp. 929–932, June 2014.[6] J. W. Kwak, H. B. Yilmaz, N. Farsad, C. Chae, and A. Goldsmith, “Two-Way Molecular Communications,”
IEEE Trans.Commun. , pp. 1–1, 2020.[7] Y. Deng, A. Noel, M. Elkashlan, A. Nallanathan, and K. C. Cheung, “Modeling and Simulation of MolecularCommunication Systems With a Reversible Adsorption Receiver,”
IEEE Trans. Mol. Biol. Multi-Scale Commun. , vol. 1,no. 4, pp. 347–362, December 2015.[8] Y. Deng, A. Noel, W. Guo, A. Nallanathan, and M. Elkashlan, “Analyzing Large-Scale Multiuser Molecular Communicationvia 3-D Stochastic Geometry,”
IEEE Trans. Mol. Biol. Multi-Scale Commun. , vol. 3, no. 2, pp. 118–133, Jun. 2017.[9] A. Noel, K. C. Cheung, and R. Schober, “Improving Receiver Performance of Diffusive Molecular Communication WithEnzymes,”
IEEE Trans. Nanobiosci. , vol. 13, no. 1, pp. 31–43, March 2014.[10] B. Li, W. Guo, X. Wang, Y. Deng, Y. Lan, C. Zhao, and A. Nallanathan, “CSI-Independent Non-Linear Signal Detectionin Molecular Communications,”
IEEE Trans. Signal Process. , vol. 68, pp. 97–112, Dec. 2019.[11] S. Giannoukos, D. T. McGuiness, A. Marshall, J. Smith, and S. Taylor, “A chemical alphabet for macromolecularcommunications,”
Anal Chem , vol. 90, no. 12, pp. 7739–7746, May 2018.[12] D. T. McGuiness, S. Giannoukos, A. Marshall, and S. Taylor, “Experimental Results on the Open-Air Transmission ofMacro-Molecular Communication Using Membrane Inlet Mass Spectrometry,”
IEEE Commun. Lett. , vol. 22, no. 12, pp.2567–2570, Dec 2018.[13] N. Farsad, W. Guo, and A. W. Eckford, “Tabletop molecular communication: Text messages through chemical signals,”
PLoS One , vol. 8, no. 12, p. e82935, December 2013.[14] L. Grebenstein, J. Kirchner, W. Wicke, A. Ahmadzadeh, V. Jamali, G. Fischer, R. Weigel, A. Burkovski, and R. Schober,“A Molecular Communication Testbed Based on Proton Pumping Bacteria: Methods and Data,”
IEEE Trans. Mol. Biol.Multi-Scale Commun. , vol. 5, no. 1, pp. 56–62, Oct 2019.[15] I. F. Akyildiz, F. Brunetti, and C. Bl´azquez, “Nanonetworks: A New Communication Paradigm,”
Comput. Networks , vol. 52,no. 12, pp. 2260–2279, Apr. 2008.[16] S. Andreescu and O. A. Sadik, “Trends and Challenges in Biochemical Sensors for Clinical and Environmental Monitoring,”
Pure Appl Chem , vol. 76, no. 4, pp. 861–878, Jan. 2004.[17] B. Alberts, D. Bray, K. Hopkins, A. Johnson, J. Lewis, M. Raff, K. Roberts, , and P. Walter, “Essential Cell Biology,” (3ed.). Garl. Press New York. , 2009.[18] B. Wang and M. Buck, “Rapid Engineering of Versatile Molecular Logic Gates using Heterologous Genetic TranscriptionalModules,”
Chem. Commun. , vol. 50, no. 79, pp. 11 642–11 644, Jul. 2014.[19] B. Wang, R. I. Kitney, N. Joly, and M. Buck, “Engineering modular and orthogonal genetic logic gates for robust digital-likesynthetic biology,”
Nat. Commun. , vol. 2, no. 508, pp. 1–9, Oct. 2011.[20] L. J. Kahl and D. Endy, “A Survey of Enabling Technologies in Synthetic Biology,”
J. Biol. Eng. , vol. 7, no. 1, p. 13,May 2013.[21] Y. Xiang, N. Dalchau, and B. Wang, “Scaling Up Genetic Circuit Design for Cellular Computing: Advances and Prospects,”
Natural computing , vol. 17, no. 4, pp. 833–853, Oct. 2018.[22] E. Bernard and B. Wang, “Synthetic Cell-based Sensors with Programmed Selectivity and Sensitivity,” in
Biosensors andBiodetection . Springer, Mar. 2017, pp. 349–363.[23] A. Tamsir, J. J. Tabor, and C. A. Voigt, “Robust Multicellular Computing using Genetically Encoded NOR Gates andChemical Wires,”
Nature , vol. 469, no. 7329, pp. 212–215, Jan. 2011.[24] T. Ellis, X. Wang, and J. J. Collins, “Diversity-based, Model-guided Construction of Synthetic Gene Networks withPredicted Functions,”
Nat Biotechnol , vol. 27, no. 5, pp. 465–471, May 2009. [25] A. Uri, An Introduction to Systems Biology: Design Principles of Biological Circuits . London, UK: Chapman & Hall,2006.[26] Y. Deng, M. Pierobon, and A. Nallanathan, “A Microfluidic Feed Forward Loop Pulse Generator for MolecularCommunication,” in
Proc. IEEE GLOBECOM , Dec 2017, pp. 1–7.[27] D. Bi, Y. Deng, M. Pierobon, and A. Nallanathan, “Chemical Reactions-Based Microfluidic Transmitter and Receiver forMolecular Communication,” arXiv preprint arXiv:1908.03441 , Aug. 2019.[28] D. Bi and Y. Deng, “Digital Signal Processing for Molecular Communication via Lego-Like Chemical Reactions-BasedMicrofluidic Circuits,” submitted to IEEE Commun. Mag. [29] G. M. Whitesides, “The Origins and the Future of Microfluidics,”
Nature , vol. 442, no. 7101, pp. 368–373, July 2006.[30] H. Bruus,
Theoretical microfluidics . London, UK: Oxford Univ. Press, 2008.[31] K. W. Oh, K. Lee, B. Ahn, and E. P. Furlani, “Design of pressure-driven microfluidic networks using electric circuitanalogy,”
Lab. Chip , vol. 12, no. 3, pp. 515–545, Nov. 2012.[32] A. G. Toh, Z. Wang, C. Yang, and N.-T. Nguyen, “Engineering Microfluidic Concentration Gradient Generators forBiological Applications,”
Microfluid Nanofluid , vol. 16, no. 1-2, pp. 1–18, Jul. 2014.[33] W. Wicke, T. Schwering, A. Ahmadzadeh, V. Jamali, A. Noel, and R. Schober, “Modeling Duct Flow for MolecularCommunication,” in
Proc. IEEE GLOBECOM , Dec 2018, pp. 206–212.[34] A. O. Bicen and I. F. Akyildiz, “End-to-End Propagation Noise and Memory Analysis for Molecular Communication overMicrofluidic Channels,”
IEEE Trans. Commun. , vol. 62, no. 7, pp. 2432–2443, July 2014.[35] R. Chang,
Physical Chemistry for the Biosciences . Herndon, VA, USA: University Science Books, 2005.[36] D. Scalise and R. Schulman, “Designing Modular Reaction-Diffusion Programs for Complex Pattern Formation,”
Technology , vol. 2, no. 01, pp. 55–66, Mar. 2014.[37] D. Harris and S. Harris,