Salinity-Based Molecular Communication in Microfluidic Channels

In this work, we present a novel portable, flexible and easy-to-use experimental setup for investigating salinity-based information transmission in microfluidic channels. At the receiver, the different salinity-levels are detected by a customized electronic circuit, which measures electrical conductivity via electrodes in the microfluidic channel. We provide a detailed description of the setup, including the microfluidic chip fabrication. Moreover, we develop a rigorous mathematical model of each testbed component and an end-to-end model of the system, which we have verified through experiments. Finally, we analyzed the error performance of the setup using optimum and sub-optimum detection algorithms, such as the Viterbi algorithm and threshold detection.


I. INTRODUCTION
M OLECULAR communication (MC) is a novel discipline in communication engineering relying on chemical molecules as carriers of information [1]. This kind of communication can be found in living organisms across all scales. For example, on a cellular scale bacteria use quorum-sensing in order to coordinate their behavior [2], [3]. In human and animal bodies hormones transmit information between different organs [4]. On a larger scale plants and animals use pheromones to communicate over the air [5]. The MC paradigm can be used to analyze living systems to provide new insights into biology, medicine and environmental sciences [6]. Besides explaining the behavior of existing biological systems, MC research also focus on the design of synthetic MC systems. This is especially beneficial in domains, where communication using electromagnetic waves is infeasible or detrimental, such as in tunnels and pipes on the macro-scale, or in the human body on the micro/nano-scale. Some of the envisioned applications for synthetic MC systems are targeted drug delivery [7] and the Internet of Bio-Nano Things (IoBNT) [8]. Even though the principle of MC has been present in nature for millennia, MC is a young research field where the first paper was published in 2005 [9]. Many theoretical works have been published in recent years, which investigate various communication and information-theoretic aspects of MC systems. While these works are very valuable for the calculation of theoretical boundaries, they typically focus on closed-form analytical solutions, which in most cases require simplifying assumptions on all components of a MC system, namely the transmitter (TX), channel and receiver (RX) [10]. Most commonly a point or volume TX are considered, both assuming instantaneous generation and release of molecules, which is hardly possible from a physical point of view [11]. Furthermore, the (advection)diffusion equation, which is mainly used to model MC channels, is a partial differential equation and thus its solution depends heavily on boundary and initial conditions. Analytical solutions can typically only be found for very simple geometries, while practical systems might exhibit more complex shapes. Finally, the RX is commonly modeled as either passive (no interaction with signaling molecules), absorbing (full interaction with signaling molecules) or reactive (chemical reaction), where the last RX type is the most practical but also most complex one [10]. In order to show that these simplified theoretical models still reflect the nature of practical MC systems, they are typically verified through simulation and experiments. Simulations are an essential part of MC research and many problems were tackled using different simulation approaches [10]. However, the limiting factor of simulations is that their outcomes rely on the underlying model assumptions, such as propagation mechanism and chemical interactions. Hence, if the importance of different effects was underestimated, simulations may still not reflect the physical reality. Experiments on the other hand do not rely on any modeling assumption and thus are an important tool for verifying and refining the assumptions made in theory and simulation. In order for experiments to be representative and reproducible the aim is to create a stable environment that reflects the physical system to be investigated. To achieve this goal various MC testbeds have been proposed in recent years. In [12] the authors proposed a tabletop spray-testbed, which for the first time demonstrated molecular communication over the air using alcohol molecules. In the This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ following years, many different types of testbeds have been proposed, which can be associated with different disciplines in science and technology. The testbeds can be classified as macro-scale (e.g., [13]) and micro-scale (e.g., [14]) and at the micro-scale a distinction can be made between biological (e.g., [15]) and non-biological (e.g., [16]) setups. A comprehensive overview of MC testbeds can be found in [17], [18]. To properly classify the testbed proposed in this work we will discuss the most related works in more detail in the following [14], [19], [20]. These works consider micro-scale non-biological testbeds and thus the proposed testbed also falls into this category. In [14], a microfluidic testbed using single-stranded DNA molecules and an integrated graphene-field-effect-transistors for the detection was presented. This testbed can be scaled down to the micro/nano-scale, which is especially important for realizing IoBNT applications. However, only low data rates (approx. 0.17 bit/s) can be achieved and the setup suffers from cross-sensitivity to non-information molecules. Due to the integration of graphene-field-effect-transistors into the microfluidic chip a complex and costly chip fabrication based on photolitography is required. In [19], a microfluidic testbed using fluorescent microbeads as information carrier was presented. The detection was performed through an optical measurements system, which does not disturb the flow of the microbeads. Compared to [14] a simple chip fabrication method based on laser ablation was used. In [20], a salinity-based information transmission through a tube network was presented. For the detection a commercial available electrical conductivity (EC) meter is employed. They identify the main characteristics of the MC channel (non-causal intersymbol interference, long delay spread, and signal-dependent noise) and develop a detection algorithm that takes these issues into account.
In this work, we propose a novel testbed for salinitybased communication in microfluidic channels. Similar to [20], we employ conductivity measurements for the detection, but we developed a customized electronic circuit that allows to be optimally adapted to the testbed and to develop a mathematical model. In contrast to the related works discussed earlier, we provide a rigorous mathematical model of each testbed component and an end-to-end system model that provides in-depth insights into the system behavior. Moreover, the model enables optimization and provides design guidelines for the testbed. Similar to [19], we employ a simple fabrication method based on 3D printing, but as chip material we use Polydimetylsiloxane (PDMS) instead of Polymethylmethacrylat (PMMA). It is important to note that due to the simple chip fabrication the testbed can be easily extended to more elaborate network of channels and, thus, provides more flexibility than the tube network presented in [20]. The proposed testbed also allows for multi-level signaling which is an important difference to the related works. We demonstrate that with 4-level signaling transmission rates up to 2 bit/s are possible. Finally, besides information transmission the testbed can be used for the precise concentration control in microfluidic chips which is important in many lab-on-a-chip applications. The main contributions of this work can be summarized as follows: 1 1 Parts of this work have been published in [21].
• We develop a novel portable, flexible and easy-to-use microfluidic testbed for salinity-based communication that can be fabricated with simple methods and easily accessible components. • We provide a rigorous mathematical model of the testbed and verify it through experimental results. • We experimentally evaluate the error performance of the testbed using optimum and sub-optimum detection algorithms. The rest of this work is organized as follows: In Section II we introduce the individual components of the proposed experimental setup. In Section III we provide a detailed mathematical model of the entire setup and in Section IV we verify the model through experiments. In Section V we evaluate the error performance of the setup and in Section VI we provide some concluding remarks and future research directions.

II. EXPERIMENTAL SETUP
In this section, we describe the main components of the experimental setup. Microfluidic pumps form the transmitter of our communication system. They apply different pressures to the microfluidic channel, which leads to different flow rates into the inlets. The two inlet flows carry saltwater and desalinated water, which mix across the microfluidic chip. Depending on the ratio of these two flow rates different salinity concentrations can be measured at the receiver. Thus, information can be transmitted by adjusting the pump pressures, which is recovered by a salinity measurement. Therefore a pair of electrodes is placed in the channel, which is connected to a conductivity measurement circuit and an analog-digital converter (ADC), which is connected to a Raspberry Pi (RPI). The RPI reads the ADC output and performs different detection algorithms presented in Section VI. A picture of the overall setup is shown in Fig. 1. The following subsections specify the most important components of the setup in more detail.

A. Pumps and Pump Control
We used the "mp6-liq" microfluidic pumps from Bartels Microelectronics GmbH [22]. These pumps are small in size (approx. 30 × 15 × 3.8 mm) and light (approx. 2g) and thus are perfectly suited for realizing a portable microfluidic testbed. According to the datasheet [22], the output pressure of these pumps depends linearly on the amplitude of the AC-voltage applied to the electrodes of the pumps. This ACvoltage is generated by the so-called mp-Highdriver driver circuit [23] also from Bartels Microelectronics GmbH. The driver circuit receives its instructions from a microcontroller via I2C and sets the pump voltages accordingly.
B. Microfluidic Chip 1) Geometry: The geometry of the channel is shown in Fig. 2. There are two inlets, which lead to a Y-junction. After the Y-junction follows a straight channel followed by a zigzag mixing structure leading to the outlet. On its way to the outlet the fluid passes the electrodes, where the salinity is measured.
The numerical values and definitions of all variables in Fig. 2 can be found in Tab. I. The variables w and h are width and  height of the channel, respectively. The variable A is the channel cross-section, L the length of the central channel (i.e., channel from the Y-junction to the outlet) and L el the distance (along the channel) from the Y-junction to the electrodes location.
2) Fabrication: The previously described channel geometry was incorporated into a microfluidic chip with the following fabrication steps (see Fig. 3).
(a) The channel structure was designed in CAD software and then 3D-printed using Polylactide Acid (PLA) as printing material. Additionally to the channel structure a surrounding frame was printed, in order to hold the PDMS poured over the work piece in the next step. The height of this frame represents the height of the final chip and was about 8mm. (b) SYLGARD TM 184 Silicone was prepared in a 1 to 10 mix-ratio between and de-aired for approximately half an hour. After pouring it over the 3D print, the entire structure rested for 24 hours for the PDMS to dry and harden. (c) A utility knife was used to carefully remove the dry PDMS from the 3D-printed structure. (d) A substrate holder of an approximate height of 2mm was printed and again covered in PDMS. This structure  was left for drying for approximately 8 hours at room temperature, such that the PDMS on the substrate holder solidifies only partly during this time period. (e) The inlet and outlet holes were drilled (diameter 1mm) and the surface of the chip was cleaned using Isopropanol. (f) The chip is placed on the substrate. The PDMS on the substrate holder is at this point just as solid, that the chip does not sink due to its own weight. Applying a slight pressure from above enables bonding of chip and substrate without the PDMS from the substrate entering the channel of the chip. The structure is then left for drying for 24 hours. (g) The microfluidic chip is cut out of the substrate holder and the electrodes are placed from above. Finally, the microfluidic chip was connected to the elastic tubes at the inlets and outlets.
C. Electrodes and Readout Circuit 1) Electrodes: We used acupuncture needles made of copper as electrodes. A copper wire was soldered to the needles and they were inserted in the PDMS at distances L el to the inlet (see Fig. 2 and Tab. I). The distance between the needles was approximately 3mm and they were placed symmetrically to the point marked as "Electrode Location" in Fig. 2.
2) Measurement Circuit: The measurement circuit serves as an interface between the microfluidic channel and the external environment. It utilizes the fact that the conductivity of water depends linearly on its salinity. We built it from discrete easily accessible components, which were soldered to a veroboard. A detailed description of the circuit is provided in Section III-C.

III. MATHEMATICAL DESCRIPTION OF THE SETUP
We start the mathematical description at the transmitter, i.e., the pumps, followed by the microfluidic channel and finally the receiver, i.e., the measurement circuit. All these components have been qualitatively described in the previous section. In this section, we present an in-depth mathematical model, which will be verified and evaluated in the next section.

A. Transmitter
We start our derivations of the mathematical model of the system by modeling the combination of pumps and highdrivercircuit by use of the information provided in the corresponding datasheets [22], [23]. The pressure drop P across the pump is related to the volumetric flowrate Q through the pump by 2 where P pump is the open-circuit-pressure drop across the pump, when its outlet is plugged and thus no flow occurs. The internal hydraulic resistance R int is a measure for the decrease of the pressure drop over the pump, if it supplies a given microfluidic setup with a flow Q. According to the datasheets, for Q = 0 the pressure drop across the pump P = P pump (see (1)) depends linearly on the amplitude of the applied sinusoidal voltage 3 V where k 1 denotes the constant amplification factor from the applied voltage amplitude to the resulting pressure drop across the pump. 4

B. Channel
Before we go into detail on the mathematical model of the channel, we qualitatively describe the working principle of the setup: The pumps inject flows of salt water and desalinated water into the respective inlets. When we assume non-negative flows (which will become important in the following) we can assume that each inlet channel, i.e., the channel from inlet to Yjunction, is filled with fluid of one of the two types. Changing the flow rates will allow injecting different concentrations of salt into the fluid streaming from the Y-junction to the outlet. Analyzing this effect mathematically will be the objective of the following discussion. In particular, we derive a pump control law, necessary conditions for the system to work properly and an expression for the steady-state salinity. Finally, we deal with the transient behavior, mainly originating from advection diffusion effects in the microfluidic channel.

1) Electrical Analogy of the Microfluidic Chip:
The behavior of fluidic systems can be described by the Navier-Stokes equations. Due to the complexity of these equations, analytical solutions can only be found for special combinations of boundary conditions and fluid properties. Thus, it is common practice in fluid mechanics to define dimensionless numbers, which allow generalization from one problem to a similar one. One important dimensionless number is the Reynolds number [25] which is a measure for the ratio of inertia to viscose drag. Thereby ρ is the fluids density, v the velocity of the fluid, a the characteristic length of the considered channel and η the dynamic viscosity of the fluid. The Reynolds number is used to determine the flow regime for a given system. In this work, we deal with laminar flow, which is defined by a low Reynolds number, i.e., typically Re < Re crit with Re crit = 2000. We will show in Section IV that this condition is met for the presented testbed. We further assume that the fluid is incompressible. Thus we can use an existing solution of the Navier-Stokes equations, that is the solution for a channel filled with incompressible fluid flowing in a laminar manner.
The following statements on the electrical analogy are based on [25]. For a circular tube the flow velocity v(r) as a function of the radial position r in the tube is given by with d the diameter of the tube, η the dynamic viscosity of the fluid, L the tube length and ΔP the pressure drop across the tube. In order to obtain a 1D model of the system, we integrate the velocity over the channel cross-section to obtain the overall flow As we will use both hydraulic and electrical resistances throughout this work, we distinguish them by including a superscript "h" for hydraulic resistances. The hydraulic resistance of a channel with circular cross-section (for example the Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. resistance of the tubes connecting the pumps to the chip) can be expressed as If the channel shape is rectangular, the hydraulic resistance is given by which is used to describe the hydraulic resistance of the considered microfluidic chip, which has a rectangular crosssection. Besides the analogy to Ohm's law (5) also an analogy to Kirchhoff's laws exists, that is the law of mass-conservation m i=0 Q i = 0 and the law of energy-conservation n i=0 P i = 0 with m and n being the number of flows Q i into a node and the number of pressure drops P i in a mesh, respectively. Based on these results, we can now proceed to apply the hydraulic analogy to our setup and analyze it using resistances and voltage sources in the way known from electrical engineering.
We consider the microfluidic channel shown in Fig. 2 and the corresponding electrical analogy in Fig. 4. In the simplified electrical circuit, the resistances R h 1 and R h 2 correspond to the input channel and the central channel, respectively. The resistance R h 1 is formed by the internal pump-resistance R h int , the tube resistance R h t,in , and the inlet channel resistance R h inlet and can be expressed as Since the internal resistance of the pump is typically significantly higher than the one of the chips inlet and that of the connecting tubes, R h 1 can be approximated by the internal pump-resistance R h int . To calculate the resistance of the central channel R h 2 , we use (7) and L, w and h from Tab. I. However, we need to account for the additional resistance from the sharp bends in the zigzag structure. According to [26], each bend contributes an additional pressure drop, for which we can account by an additional term that we add N times for the N bends of the zigzag. The pressure drop across one bend can be calculated as [26] with ζ = 1.3 for a 90 • bend and b ≈ 668 Re −1 , allowing the calculation of a single bend as By neglecting the resistance of the outlet-tubing R h t,out , R h 2 can be calculated as Note that neglecting the outlet tubing (diameter 1mm, length approx. 20cm) is justified, since for the given geometries its contribution is small (see (6)). We now assume that the lower inlet-channel (see Fig. 2) is entirely filled with salt-water (salinity s 1 > 0) flowing with a velocity v 1 = Q 1 /A and the upper inlet channel is entirely filled with desalinated water (salinity s 2 = 0) flowing with velocity v 2 = Q 2 /A. At the junction the two fluid streams meet, but do not mix. While propagating across the channel and aided by the mixing structure, they merge into each other. At the receiver we assume the fluids are well mixed.
If both input flows are held constant for a sufficiently long time, the fluids salinity s el will have converged to a static value. This value can be simply determined by the massbalance equation [27]. This equation states that the amount of salt carried by the two inlet streams needs to be also contained in the outlet flow. This means that Q 1 s 1 = s el (Q 1 +Q 2 ) and consequently In order to have a linear dependence of the salinity at the electrodes on the flow induced by the pumps, we demand that It can be shown using the electrical analogy that this requires the sum of pump-pressures to be constant, i.e., which is obviously satisfied by choosing and Thereby, P ref = k 1 V 0 denotes the chosen offset, which is constant over time. It is necessary, since the pump pressures can only take on positive values. The voltage V 0 is the pump input voltage required to realize the open circuit pump pressure P ref and ΔV i is the voltage used to encode the i-th symbol.
The hydraulic analogy allows the calculation of the flows Q 1 and Q 2 as follows (see Fig. 4) and Substituting (15) and (16) into (17) and (18) we obtain At the beginning of this section we highlighted that Q 1 and Q 2 are not allowed to get negative. The reason for this is that, if one of the flows would become negative, we would pump salt-water into the reservoir for desalinated water or vice versa, which would obviously disturb the entire working principle of the system. From (19) and (20), however, we notice that, even if we apply non-negative pressures P 1 and P 2 , one of the flows can still become negative, if This allows to derive a maximum allowed voltage difference at the pumps as This is the available voltage range in which we can encode information symbols, e.g., for a four symbol alphabet this range is divided into four levels. Moreover recall that R h 1 and R h 2 are mainly defined by the pump and the central channel resistance, respectively. The pump resistance is constant, while the hydraulic resistance of a annular tube (see. (6)) or quadratic channel (a special case of (7)) scales approximately with the inverse fourth order of the channel diameter (for a rectangular channel the formula (7) is a bit more complex, yet shows a similar tendency). Consequently for large channels, the entire input range is usable (R h 1 R h 2 ), while shrinking the channel leads to a reduction in the usable input range and thus also a reduction in the average signal energy. We will revisit this idea after dealing with the transient behavior in the following.
2) Advection-Diffusion Channel: In the previous subsection we assumed that the two fluids mix across the channel. For a straight channel with laminar flow, however, this will only take place, if the Brownian motion happening within the two fluids has sufficient time to carry particles radially across the channel. The time the fluid resides in the channel, however, depends on the flow velocity of the fluid. As already mentioned prior, it is common practice in fluid mechanics to distinguish scenarios by dimensionless numbers. The number describing the ratio between residence time in the channel and time taken to diffuse across the channel is the Péclet number [25] with v the longitudinal flow-velocity, D the diffusion coefficient and a the characteristic length of the channel. Obviously, if flow dominates over diffusion, then Pe 1 and if Pe 1 diffusive phenomena dominate the process. Thus, efficient mixing requires a low Péclet-number. We will show in Section IV that in absence of a mixing structure the Péclet number will be much larger than 1. This has the consequence that without additional mixing-structures we cannot expect mixing to happen. Many mixing structures for microfluidic channels have been proposed [28]. We decided for the zigzag structure proposed in [24], due to the simplicity of integrating it on the chip and due to its passive working principle. The mixing structure can be regarded as a device that increases the effective diffusion coefficient of the entire structure. This can be seen from the fact that at the outlet of a mixing channel the fluid is mixed as well as it would have been, if diffusion had sufficient time to mix the fluids. Consequently, we cannot neglect diffusive effects, if we have a mixing structure on our chip. In other words, the presence of the mixing channel artificially decreases the effective Péclet number. Therefore, we use the advection-diffusion equation to model the transient behavior of the system ∂s(r , t) with r a spatial vector, t the time, D the diffusion coefficient, v the velocity vector and s(r , t) the salinity as a function of space and time, respectively. The first term on the righthand side reflects diffusive effects, while the second term reflects advective effects. Based on the previous discussion, we can view the mixer as a structure that reduces the axial distance to be traveled to reach the same degree of mixing as in a much longer straight channel. For a sufficiently long tube (or a shorter tube including a mixing structure), however, the 3D model (24) can be accurately approximated by a 1D model [29]. This phenomenon is known as Taylor-Aris dispersion. Its most fundamental finding is that after a sufficiently long distance traveled in a tube, the fluid concentration profile will always tend towards a bulk flow profile and the crosssection averaged concentration (or salinity) in response to an impulsive concentration burst at the origin is given by The Taylor-Aris dispersion model predicts that this block will spread axially a lot faster than one would suspect from Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. the regular diffusion coefficient D. In fact, according to this theory, where a is a characteristic length and K denotes a constant describing the geometry. The impulse response h ad (t) of this system evaluated at a distance z = L el from the origin is then given by Note that this is just the advection-diffusion channel impulse response of the microfluidic channel and does neither include the encoding process (discussed in Section III-B1) to get from an applied voltage to a salinity, nor the conversion of the salinity to the voltage at the receiver (see Section III-C).

3) Concluding Remarks:
In this section, we derived the model of the microfluidic channel. In the following, we will discuss some interesting properties of the model. From (22) we concluded that the maximum usable input voltage range and thus maximum average symbol energy is obtained, when the channel diameter is large. However, (26) suggests that the effective diffusion coefficient increases with increasing channel diameter. The effective diffusion coefficient is a measure for the spread of the impulse response (44) of the system and thus for the strength of intersymbol interference (ISI) for a given symbol duration. Consequently, there is a trade-off between signal energy and higher intersymbol interference. Based on the mathematical description several degrees of freedom exist (e.g., channel diameters, length of the inlet and outlet channels, pumps etc.). In this work we chose the geometry to observe all previously described phenomena, such as limited allowed input voltage, intersymbol interference etc..

C. Receiver
The receiver has the task of measuring and interpreting the information encoded in the salinity at the electrodes. Therefore, the circuit shown in Fig. 5 was designed, which will be discussed in the following.
We know from [30] that the conductivity σ(s) of salt-water depends on its salinity s according to with a constant conversion factor α. In a practical setup, however, the presence of metallic electrodes leads to the occurrence of additional phenomena that might influence the measurement [31]. In particular, when there is an interface between a metallic electrode and an electrolyte solution, a so-called double layer will form. This double layer can be modeled by a capacity C D that is in series to the resistance of the fluid R D (this is a form of the so-called Randles circuit, shown in Fig. 6).
Consequently, DC-current is blocked from flowing through the element, which makes determination of the conductance in (28) impossible. Therefore, the proposed measurement circuit (see Fig. 5) uses an AC-measurement principle. The device-under-test (DUT) is formed by the element shown in Fig. 6 and has the impedance If we apply a sinusoidal voltage u D (t) =V D sin(ω D t) to this system and choose ω D sufficiently high, so that R D j ω D C D , the current through the DUT can approximated through leading to the output voltage of the operational amplifier given by where g D (s(t)) = 1/R D (t) = k 2 s(t) with k 2 a constant depending on the geometry and location of the electrodes. The salinity s(t) will depend on the time. Since changes of s(t) are a consequence of advection-diffusion effects, these changes can be assumed to happen slowly compared to the changes of the sinusoidal voltage applied to the DUT. We can thus consider (31) as a baseband signal amplitude-modulated onto a sinusoidal carrier of angular frequency ω D . We rectify this signal, apply the Fourier transform and use the knowledge that g D (t) ≥ 0 to obtain with the convolution-operator * and G D (ω) the Fourier transform of g D (t).
The second part on the right had side of (33) has a line spectrum, which can be calculated using the Fourier series as Thus, the convolution in (33) yields .
Since we assume that G D (ω) only contains spectral components in the range from ω = 0 to ω = ω S ω D , we can recover it by passing it through a low-pass-filter with a corner frequency ω C with ω S ω C ω D . Ideal low-pass filtering would give at the receiver output. The corresponding time-domain signal is given by with k 3 = k 2 RV D 2 π . Note that we have two ways to alter k 3 , namely by changing eitherV D or R. Due to the fact that Randles model only holds for sufficiently small voltages, we keepV D constant and set the measurement circuit gain by changing R.

D. End-to-End Model
In the previous subsections, we dealt with the physical motivations and derivations governing the setup. Here, we connect all these derivations in order to obtain an end-to-end model. We again start at the transmitter side and consider a symbol u i , which is mapped to a voltage ΔV i . A message x(t) consisting of N symbols can then be written as with As mentioned previously, we assume that the pumps instantaneously encode the concentration s el at the channel input (even though at the input the concentration is actually encoded in two parallel phases, which are not mixed). Then, according to (12), each of the individual rectangular pulses in (38) leads to an instantaneous encoding of a salinity s i = s 1 Q 1 Q 3 . Substituting (19) and (20) into this formula gives Note that the first expression is always the same, no matter which symbol is encoded. This expression reflects the constant offset imposed by V 0 . The second term depends on the actual symbol, since it depends on ΔV i . Since the microfluidic channel output can be obtained by convolution of x(t) and h ad (t), we get At the electrodes the signal y ad (t) is then measured and converted into the electrical signal y(t), which can be written as Introducing the constants k 4 = k 3 s 1 2 and k 5 = k 3 s 1 B 1 2A 1 , we obtain the end-to-end model of the entire system as

IV. EXPERIMENTAL RESULTS
In this section, we verify the theoretical models presented in Section III through experiments. 5 Due to limitations in space, we will focus on the most important parameters, which can be distinguished into three categories: A. Parameters that are known from the literature or were designed by us, for example η, L el . B. Parameters that are unknown but can be directly or indirectly measured with the available equipment, for example v, R h 2 . C. Parameters that cannot be measured but estimated from a specific measurement, for example D eff . In the following, we start by providing the most important numbers from the literature, our chip design and datasheets. Then, we introduce the used measurement equipment and explain the different measurements taken. Finally, we deal with those quantities that cannot be measured but can be estimated from signals obtained by measurements. We conclude the section by showing that the deterministic model accurately describes the average system behavior and by analyzing the predominant noise source.

A. Parameters From the Literature 1) Geometrical Constants:
In order to calculate the hydraulic resistances and for the behavior of the advectiondiffusion channel the geometries of the channel are required, which are provided in Tab. I. However, we didn't specify the characteristic length a and the constant K. The characteristic length a for a rectangular duct is given by [14], [32] 5 It is important to note that in the following we verify the models for the parameters of a specific setup. However, during our investigations we have used different setups for all of which we observed a good agreement with the proposed models. The constant K for the Taylor-Aris dispersion (see (26)) is approximated by the parallel plate scenario. We obtain from [33] 2) Fluid Properties: Tab. II summarizes the fluid properties, most of which can be found in the literature. All values vary slightly with salinity and temperature. The values in the table were therefor chosen within the possible range of variation 3) Pump Resistance: As discussed in Section III, the hydraulic resistance R h 1 is dominated by the internal resistance of the pumps R h int , which can be obtained from the datasheet [22]. In analogy to electrical engineering, the internal resistance is the ratio of open circuit pressure P oc to shortcircuit flow Q sc , which are both given in the datasheet. We choose a frequency of 100Hz as the excitation frequency of the pumps provided by the highdriver circuit. For this frequency we obtain P oc = 500 mBar = 5 × 10 4 Pa and Q sc = 0.005 L/min ≈ 8.3 × 10 −7 m 3 s −1 , resulting in The factor k 1 , which is also needed to describe the pumps behaviour, can be modeled by with ΔV oc a voltage step and ΔP oc the corresponding pressure step, which are obtained from the datasheet of the pump [22], [23].

B. Measured Parameters
Next, we characterize the hydraulic part of our mathematical model. Note that for R h int we did not care about the flow regime we are in, since the datasheet provided the internal resistance. However, in order to calculate R h 2 , we need to verify our assumption that we are in a laminar flow regime to ensure that the presented formulas can be applied. This means we need to determine the Reynolds and Péclet number.

1) Reynolds and Péclet Number:
To derive these numbers the only unknown quantity is the flow velocity v. To determine it we use the liquid flow sensor and the accompanying software from Sensirion [37]. Since we know that the flow rate from the Y-junction to the outlet is symbol-independent, we chose to measure the flow rate when both pump voltages are set to V 0 = 125 V, which lies in the center of the usable input voltage range from 0v to 250v. From (19) and (20) we notice that the flow we measure at the outlet is just Q 3 = 2A = 360 µL s −1 . Therefore, we obtain a flow-velocity v = 2A wh = 0.0323 m s −1 yielding a Reynolds number of which means that our flow-regime is laminar. Knowing the velocity also enables calculating the Péclet number Pe hyp We use the term "hyp" to denote that this is the Péclet number we would hypothetically get without a mixing structure. Consequently, without a mixing structure, no mixing would occur.
2) Hydraulic Resistance from Y-Junction to Outlet: Having shown that the assumption of laminar flow is valid, we compare the theoretical value for the resistance R h 2 (see (11)) to the measured one. This requires a cut through the chip at the Y-junction, 6 to isolate the resistance R h 2 . 1) For this experiment we used microfluidic pumps [38] from Fluigent. These pumps are more suitable for this measurement, since they have an integrated pressure control. Furthermore, they provide a measured outputpressure, which will be needed for the further calculations. These pumps work differently than the ones we used previously. They supply pressurized air to a liquidfilled reservoir, pushing the fluid out of the reservoir.
In the steady-state, the equivalent internal resistance is only determined by the resistance of the tubing from this reservoir to the chip, which is not negligible in this case, since this tubing has a significantly lower diameter than the one previously used. 2) At the output of the pump reservoirs we attached the Sensirion flow sensor. Since it has a small hydraulic resistance compared to R h 2 , we neglect its influence on the measurement. 3) We choose a pressure P test = 100 mBar and measure the flow Q test,0 through the tubing with no chip attached. 4) Using a scalpel we separate the chip at the Y-junction and use the outlet of the Y-junction as inlet for this experiment. Then, we connect it in the same manner as described above. 5) We apply the same pressure P test and observe how the flow changes. We call the newly measured flow Q test,1 . 6) Inspired by hydraulic analogy (see Fig. 7) we can calculate the hydraulic resistance as Please note that this procedure is carried out 10 times and the results were averaged. We obtained an estimate for the resistance R h 2,meas = 9.42 × 10 10 Pa s m −3 , which is in a good agreement to the calculated value R h 2,calc = 1.01 × 10 11 Pa s m −3 obtained from (11) using the parameters presented in the respective tables. Electrical analogy of the hydraulic measurement procedure. Measurement of closed-circuit-flow (left) and when the resistance to be measured is connected (right).
With this finding, we are now able to predict the usable input voltage. Using the previously determined values for R h 1 = R h int = 6 × 10 11 Pa s m −3 (see (47)) and R h 2 ≈ 1 × 10 11 Pa s m −3 we conclude that, according to (22), the maximum usable voltage difference is ΔV max,theoretical = 3 4 V 0 . In order to stay well within that range, we our practical voltage difference much smaller than this theoretical value. Choosing V 0 = 125 V, we use from now on ΔV max,theoretical = 37.5 V. This means, to encode a 4-level symbol, the voltage ΔV i (see (43)) can take on the values [0.0 V, 12.5 V, 25 V, 37.5 V] and in the 2-level case the values [0.0 V, 37.5 V] with each voltage corresponding to a specific symbol. As discussed in Section III-B1, each voltage can be related to a specific encoded concentration level and, thus, the encoding scheme corresponds to concentration-shiftkeying (CSK) [39].
3) Measurement Circuit: Finally, the variable k 3 (see (37)) relating the salinity in the channel to the measured voltage is also a measurable quantity. Thus, in the following, we will show an exemplary measurement of this quantity. Therefore, we prepared syringes filled with saltwater at four different salinity-levels: s 0 = 2.5 g L −1 , s 1 = 5 g L −1 , s 2 = 7.5 g L −1 and s 3 = 10 g L −1 . By use of syringes, the different liquids were injected in the order presented and then in reversed order. Fig. 8 shows the outcome of the measurement. The measurement starts with an air-filled channel. Thus the first step is higher than the others, since the model proposed does not hold for this case. Then the different fluids are injected. The step size is always approximately equal. The sharp drop between two injections reflects the process of changing the syringe, where for a short time some air enters the channel and, thus, the measured value drops. Obviously, this is the reason, why for the first and last step this drop is not visible. From the shown measurement, k 3 was deduced by dividing the change in voltage (approx. 0.12v) by the change in salinity (approx. 2.5 g L −1 ) resulting in k 3 ≈ 0.048 V L g −1 . Consequently, we have shown that for a given electrode placement it is possible to determine the constant k 3 , which also means that we can translate the measured voltage signal to an actual salinity level in the channel. This becomes especially important, if we want to use the measurement circuit as an interface between the MC environment and the surrounding world. At this point we would like to remind the reader, that for the experimental section of this work multiple chips where used. While the hydraulic and diffusive behavior does not vary to strongly between chips, the constant k 3 , which is dependent on the placing of the electrodes, changes significantly. However, as can be seen in (44) it is only a scaling factor for the salinity signal. Thus, in order to make comparing measurements on different chips easier, in the following, we will shift our measurements to have zero offset and re-scale them. Whenever we use this approach, we refer to the voltage as "normalized voltage".

C. Estimated Parameters
Besides the flow-velocity v, the effective diffusion-coefficient D eff is the most important quantity affecting the transient behavior of the system. We can, however, not measure this parameter with the available measurement instruments and therefore need to estimate it. We proceed as follows: 1) We collect multiple impulse responses. This is done by sending 100ms pulses with an amplitude of 37.5v and recording the voltage at the receiver. The recorded voltage is then again shifted and rescaled, since we are only interested in the temporal dynamics of the system and not in the gain or offset. 2) We use the MATLAB optimization toolbox to fit the average impulse response to the advection-diffusion equation (27) using the Levenberg-Marquard algorithm. 3) Since the problem under consideration is a non-convex optimization problem, in general it has multiple local solutions. To check the plausibility of our obtained value, we consider also the velocity v as a degree of freedom. Then, we compare the value obtained from the estimation to the one we measured. Following this procedure, we obtained the plot shown in Fig. 9. The flow velocity obtained from the optimization algorithm is v = 0.055 m s −1 , which is reasonably close to the measured value of v = 0.032 m s −1 . Thus, we are confident that we have found the desired solution. The algorithm estimates D eff = 1.24 × 10 −4 m 2 s −1 . Using the measured velocity v in (26) we obtain D eff = 4.4 × 10 −4 m 2 s −1 . Again, we notice a good agreement between the results. 7 We notice a  significantly higher variation in the measurement of the average impulse response than we will see in Figs. 10 and 11. This is due to the fact that we used a short (100 ms) pulse of 37.5 V to mimic a Dirac pulse. Obviously the signal energy is thus very little. However, as will be shown in Section IV-E, the noise power is amplitude independent. Thus the signal-tonoise ratio (SNR) is quite low in this plot compared to the other plots, where the energy of the input signal (and thus the output signal energy) is much larger.

D. End-to-End Model
Next, we want to show that the end-to-end model (43) actually reflects the measured signal in a typical transmission scenario. For this we consider two different scenarios: 2-level and 4-level transmission with each a symbol duration of 0.5 s. For each scenario we sent a random sequence of 2-or 4-level symbols, respectively, into the real system and measured the voltage at the system output. We compare this to the outcome of the end-to-end model, when the same input From these results we conclude that the proposed model approximates the physical system reasonably well.

E. Noise
So far, we did not consider the noise present in the system. The noise arises due to various non-ideal behaviour of all involved component. For the pumps, we have switching uncertainties, variations in the injection amount and pump oscillation-induced noise. In the channel we get noise from the temperature-dependence of the hydraulic resistances, eddy phenomena at the mixing structures and noise due to variations in the particle concentration (so-called counting noise [40]). Finally, at the receiver we get noise due to the coupling between the electrical field between the measurement electrodes and the concentration profile, as well as noise from the electronic components and ADC quantization noise. Please note that a peculiarity of the counting noise is that its amplitude is proportional to the number of molecules. In terms of concentration this means, that the measured signal would exhibit higher noise variance for higher concentrations. This is fundamentally different to most noise sources in electrical engineering, where the noise power is typically independent of the signal [41]. This leads to the question, whether the dominant form of noise in our testbed is signal-dependent or independent. Therefore, we carry out the following experiment. We send different salinity-levels and keep the symbol duration long enough, so that a steady state will occur. Then we run a moving variance filter across the collected data. If the variance for different steady states is always the same, we can conclude that the variance is independent of the signal and vice versa. Fig. 12 shows the outcome of the experiment. The first subplot shows the received signal. The moving variance filtered signal is only meaningful in those time frames where the received signal reaches a static value. Therefore, the third subplot shows a zoomed view of the moving variance filtered signal on those points in time where the received signal has a constant value. We notice that the noise is signal independent, since we do not see a significant difference in the signal variance when the received signal has different values. This result is plausible since already the thermal noise in one 10 kΩresistor (of which we have many in the measurement circuit) at room temperature and for a bandwidth of B = 1 kHz yields a noise voltage of σ r = √ 4k b TR noise B = 8.13 × 10 −7 V [41]. If we relate this value to an average signal voltage of 0.5v the ratio between average signal voltage and noise variance is 1.6 × 10 −6 . On the other hand, if we relate the molecule counting noise 8 [40] to the average concentration signal at a salinity of 10 g L −1 we get a ratio of 1.825 × 10 −8 , which is approximately two orders of magnitude smaller, than the ratio for the thermal noise produced in one single resistor. Note that we assumed a cuboid receiver with d el the longitudinal distance of the electrodes in the channel, which to some extent is a simplification, since our receiver is not a simple counting receiver, but interacts with the molecules. Yet, it can be clearly seen, that the dominant noise cannot stem from the molecule counting process in our setup. Thus we assume signal independent noise for our testbed with the following variance, which can be directly obtained from the third subplot in Fig. 12 Now that we know that the noise is signal independent, we can provide the final end-to-end model as with n(t) ∼ N(0, σ 2 n ) and y(t) from (43). Recall, that for better comparability between experiments, we use normalized versions of the signals above. In order to add the correct noise amplitude σ n,norm to the normalized signal y norm (t) we define the factor κ = Δy σn = Δynorm σn,norm , with Δ indicating the peak-to-peak amplitude. Fig. 12 demonstrates a typical operating point of the setup, for which Δy ≈ 0.4 V and the σ n ≈ 10 −3 V yielding κ = 400. Finally, Tab. III provides all necessary numerical values, to be able to simulate the system dynamics and noise behavior.

V. COMMUNICATION PERFORMANCE
In this section, we analyze the performance of the proposed testbed from a communication perspective. Due to the fact, that data-collection in MC-systems is time consuming, two approaches for the error performance evaluation exist which we refer to as "CIR-based" and "experiment-based". In the CIR-based approach, for example utilized in [42], the CIR of the system is estimated from data and used as a simulation model to evaluate the error performance. Due to the fact, that these simulations are fast compared to the real system, arbitrarily many sent symbols can be simulated and the error performance of the idealized system can be accurately determined. The experiment-based approach on the other hand relies on data measured from the real setup. Hence, in the same time only a limited number of samples can be collected and, thus, this approach provides a first idea of the error performance of the system. However, the major advantage of this approach is that captures all non-idealities that are not considered in the CIR-based method. Since we want to keep this analysis closely tailored to the practical system we go for the experiment-based approach. In the first subsection, we present the data collection and how the data is pre-processed. Then, we present the used detection algorithms and finally we demonstrate the error performance of the different detection algorithms on the measured data.

A. Data Collection and Preprocessing
We are interested in the best achievable communication performance for the proposed testbed. Since some of the used detection algorithms require channel knowledge, we need to choose a way to provide the channel impulse response (CIR). Basically, there are two possibilities: 1) Estimate the average CIR in advance (see Section IV).
2) Estimate the CIR for the actual transmitted data sequence. Each approach has its advantages and disadvantages. Using the average CIR allows detection with less delay, but worse performance, since the CIR can change slightly when the data were not all collected at once. On the other hand, estimating the CIR from the actual data introduces more delay, since we need to collect sufficient symbols to estimate a meaningful CIR, before we can start the detection process. Since our aim is to evaluate the best-case performance of the testbed, we chose the second approach. The workflow for the error performance evaluation of the testbed can be described as follows: where L is the length of the considered FIR-filter and M is the number of samples in the sequence. Then this impulse response is fitted to the CIR in (44) in another optimization step. The fitted impulse response is stored in the vector h fit . 3) Detection: Application of different detection algorithms on each dataset. 4) Performance evaluation: Computation of the number of symbol errors for each detection algorithm and dataset.

B. Detection Algorithms
In this section, we introduce various detection algorithms, namely the maximum likelihood sequence estimator (MLSE), slope detector and threshold detector. The MLSE provides the optimal performance, but requires a lot of computational resources. We use it, to provide a benchmark performance. However, in a typical MC environment (especially for in vivo applications) the computational resources are limited. Thus, for the practical application we consider slope and threshold detection as low-complexity alternatives.
1) Maximum Likelihood Sequence Estimator: 9 The detected symbol sequenceû is determined by the following detection rule [43] where p(y |ũ) denotes the likelihood of observing the received signal y for a trial symbol sequenceũ. Since the channel is characterized by additive Gaussian noise (see (54)), the decision rule can be simplified as follows [43] u = argmiñ u y − w ũ os (ũ), h fit 2 , whereũ os (ũ) results from oversamplingũ and w corresponds to the deterministic system response, which is obtained through convolution of the input signal sequencẽ u os and the channel impulse response h fit . For the practical implementation of the MLSE, we use the Viterbi algorithm presented in [20], which accounts for non-causal intersymbol interference (a property of the advection-diffusion channel) and reduces the complexity by limiting the number of states.
2) Slope Detection: The slope detection is a lowcomplexity alternative to the MLSE and is based on the increase/decrease of the received signal. We employ the following heuristic for the detection of the pth binary symbol where ε denotes the detection threshold and Δy can be calculated as follows The output at the end of the symbol period is thus compared with the value l sampling periods before. A simple explanation for utilizing this detection scheme can be given as follows.
Since our system input consists of a sum of step functions, the system output for a sequence of the same symbol remains in the same range. Therefore, if the change of the system output is smaller than ε, the same symbol is detected again. If, on the other hand, the previously detected symbol is a binary "1" ("0") and the signal decreases (increases) by more than ε, a binary "0" ("1") is detected. This detection method can be easily extended to higher order modulation in a straightforward way, but for the sake of clarity only the heuristic for binary modulation is given here. For a specific symbol duration, the parameters ε and l are determined via a grid search with a simulation model before decoding the actual data. It is important to note that for a small symbol duration, i.e., when the symbol duration becomes smaller than the time at which the impulse response reaches its maximum, this detection scheme no longer works with a reasonable performance.

3) Threshold Detection:
The threshold detector is the simplest and most intuitive detector. In the binary case, it compares the received value to a threshold. If the received signals amplitude is below the threshold the detector outputs a binary "0" otherwise a "1". In the multi-level case, multiple thresholds are defined and the detector assigns each signal amplitude range to a different symbol. This detector is the computationally most inexpensive one and, thus, very well suited for MC applications. However, it cannot deal with ISI, which limits its performance.

C. Error Performance Comparison
In the following, we present the error performance of the system. Tabs. IV and V summarize the number of symbol errors of the different algorithms for different symbol alphabets and symbol durations. As expected, the Viterbi algorithm performs best for all scenarios. For a symbol duration of 1 s all algorithms achieve zero errors for binary transmission. Moreover, for the binary case the Viterbi algorithm achieve zero errors down to a symbol duration of 0.5 s. For the 4level transmission more errors occur and for high ISI scenarios (i.e., low symbol duration) the performance of low complexity detectors degrade significantly. Finally, we notice that the Viterbi algorithm provides an error free transmission at a rate of 2 bit/s. Although we only transmitted 100 symbols this indicates that the proposed testbed provides a reasonable information rate for MC systems.

VI. CONCLUSION
In this paper, we presented a novel eas-to-use testbed for molecular communication on a microfluidic chip based on the salinity of a fluid. The testbed components of the experimental setup are cheap and easily accessible. We presented detailed mathematical models of all testbed components and connected them to the physical and geometrical properties of the system. Since the testbed can be accurately described mathematically, this can be used to design a customized testbed for experiments across different domains of molecular communication. Finally, the communication capabilities of the testbed were evaluated in terms of error performance using different detection algorithms. Future works will include the investigation of multi-molecule transmission, machine-learning-based detection and the extension of the receiver towards on-chip electrodes.