IEEETRANSACTIONSONNANOBIOSCIENCE,VOL.15,NO.7,OCTOBER2016 765 A Biochemical Oscillator Using Excitatory Molecules for Nanonetworks ∗ Ethungshan Shitiri, Student Member, IEEE, and Ho-Shin Cho, Member, IEEE Abstract—For nanonetworks to be able to achieve large-scale With the advantage of being minuscule, the applications functionality, such as to respond collectively to a trigger, syn- of nanonetworks are envisioned over a wide range of fields: chrony between nanomachines is essential. However, to facili- for quality control and to develop new materials in industrial tate synchronization, some sort of physical clocking mechanism applications, to control pollution in environment monitoring is required, such as the oscillators driven by auto-inhibitory molecules or by auto-inducing molecules. In this study, tak- applications [1], to deliver drugs in nanomedicine applica- ing inspiration from the widely studied biological oscillatory tions[2],andtodetectnuclear,biological,andchemicalagents phenomena called Calcium (Ca2+) oscillations, we undertake a in military applications [5]. However, being minuscule limits different approach to design an oscillator. Our model employs the capability of the nanomachines and would require them threedifferenttypesofexcitatorymoleculesthatworkintandem to be coordinated or synchronized in order to achieve large- to generate oscillatory phenomenon in the concentration levels scale functionality, such as releasing drugs in coherence to of the molecule of interest. The main objective of the study is to model a high frequency biochemical oscillator, along with a targeted site. In this regard, the pertinent issue here is, the investigations to identify and determine the parameters that however, not about how synchronization can be achieved, but affect theperiodoftheoscillations.Theinvestigationsentailand most importantly the issue of how to implement a physical highlightthedesignofthereserveunit,areservoirofthemolecule timing unit [6], such as a clocking device or an oscillator of interest, as a key factor in realizing a high frequency stable biochemical oscillator. within a nanomachine, so that the nanomachines can simply synchronize their oscillators in order to be synchronizedwith Index Terms—Biochemical oscillator, calcium oscillations, each other. In that way, the oscillator(s) can facilitate the molecular communication, nanonetworks, nanomachines. process of synchronization. Until now, only a couple of studies have addressed this I. INTRODUCTION issue. In [7], auto-inhibitorymoleculesare employedto make MOLECULAR nanonetworks or simply nanonetworks nanomachines oscillate. Auto-inhibitory molecules have a are defined as the internetworking of nanomachines peculiar behavior of inhibiting the chemical process that is based on, but not limited to, molecular communication [1]. responsible for their generation. This inhibiting behavior is Molecular communication employs molecules as information triggered when the auto-inhibitor’s concentration increases to carriers [2] and, over the years, it has been garnering a an upper threshold, and, it subsequently causes a decrease in lot of interest over other forms of nanoscale communica- their concentration. As the concentration decreases and falls tion techniques such as the traditional electromagnetic and below a lower threshold, the inhibition is stopped and the acoustic communication, and a not so traditional contact- generationprocessisreinitiated,allowingthe concentrationto based communication called nano-mechanical communica- increase. These processes keep repeating, forming the oscil- tion [3, p. 177]. What makes molecular communication so lations. In a more recent study [8], auto-inducing molecules attractive for nanonetworks are for the facts that it is natu- are employed to make a group of nanomachines oscillate. rally designed for the nanoscale, compatible with biological Unlike the auto-inhibitors, auto-inducers facilitate their own systems, and energy-efficient [4]. Nanomachines, which can generation.Whenthenanomachinessensetheauto-inducersin be biological or electromechanical, are the functional units the environment, they generate and release the auto-inducers of a nanonetwork that are capable of performing basic func- until the concentration of the auto-inducers reaches an upper tions of sensing, actuating, computing, and communicating threshold.Oncetheconcentrationreachestheupperthreshold, (information transfer) with one another [1] [4]. the nanomachines stop the release of the auto-inducers and Manuscript received August11,2016;revised October 14,2016;accepted this action causes a decrease in the concentrationof the auto- November 1, 2016. Date of publication October 18, 2016; date of current inducers. As the concentration falls to a lower threshold, the version December 20, 2016. This study was supported by the BK21 Plus nanomachinesbegintoreleasetheauto-inducerscausingarise project funded by the Ministry of Education, Korea (21A20131600011) Asteriskindicates corresponding author. intheconcentrationagain.Theresultantperiodicincreaseand E.ShitiriiswiththeMobileCommunicationLaboratory,SchoolofElectron- decrease in the concentration of the auto-inducers constitutes icsEngineering,KyungpookNationalUniversity,Daegu41556,SouthKorea. ∗ the oscillations. H.-S. Cho is with the Mobile Communication Laboratory, School of ElectronicsEngineering,KyungpookNationalUniversity,Daegu41556,South The basis of the first oscillator [7] is from a naturally Korea(e-mail: [email protected]). occurringphenomenon,wherein,thecouplingoftwofeedback Color versions of one or more of the figures in this paper are available signals, excitatory and delayed inhibitory, results in oscilla- onlineathttp://ieeexplore.ieee.org. Digital ObjectIdentifier 10.1109/TNB.2016.2616539 tions [9]. Similarly, the basis of the second oscillator [8] is 1536-1241©2016IEEE.Personaluseispermitted, butrepublication/redistribution requires IEEEpermission. Seehttp://www.ieee.org/publications_standards/publications/rights/index.html formoreinformation. 766 IEEETRANSACTIONSONNANOBIOSCIENCE,VOL.15,NO.7,OCTOBER2016 also from a naturally occurring phenomenon, called quorum sensing[10].Infact,severalbiologicalactivitiesexhibitoscil- latorybehaviorandarewidelystudied.Theperiodicsignaling of cAMP (cyclic adenosine monophosphate) between cells, periodicincreaseanddecreaseincalciumconcentrationinside a cell, and circadian rhythms are few examples of biological oscillations [11]. Such systems are commonly referred to as biological oscillators. Biological oscillators, in general, have been studied for decades and in recent years there have been several attempts and breakthroughs in designing biological oscillators, such as genetic oscillators [12] [13] and chemical micro-oscillators [14]. These oscillators are designed from a biologicalsystemengineerperspectiveandtendtohaveacom- Fig. 1. Working dynamics of Ca2+ oscillations found in astrocytes. Solid mon feature—large period of oscillations. Such feature may lines indicate the transport/feedback mechanisms and dashed lines indicate not be desirable for nanonetworks that are deployed in time- the fluxes. Plus sign indicates positive feedback and minus sign indicates critical applications. The oscillators proposed in [7] and [8] negativefeedback. (ThefiguredepictsonlythemajorsignalingpathsofIP3, Ca2+,andPLCδ.) are probably the first studies on biological oscillators (with comparativelysmaller periodsof oscillations) from a commu- nication engineering perspective. However, in [7] there is no signals is reduced. To validate this statement, we provide an clarity on developing the oscillator as a separate functional elaborate discussion in Section III-C. unit within the nanomachine, while in [8] there is no such The remainder of this paper is organized as follows. attemptatallas a groupof nanomachinesis consideredto act Section II provides an overview of the signaling principles as an oscillator. of Ca2+ oscillations. Section III presents the working of the Giventhestochasticpropagationnatureofthebiomolecules, proposed biochemical oscillator along with the mathematical viz. Brownian motion, it is imperative to design the oscillator framework,andtoshowthatahighfrequencyoscillatorcanbe as a separate functional unit such that the oscillations, whose obtained,theconceptofanoperatorisintroducedhereaswell. molecules are already impaired by Brownian motion, are SectionIVdiscussesthesimulationresultsandtheparameters not interfered by unwanted noises. Noises, here, refers to thataffectthe meanandvariancein the periodof oscillations. the residual or lingering molecules from other processes that Finally, Section V concludes the paper. Beyond this point, can interfere with movements of the oscillatory molecules. the proposed biochemical oscillator will be interchangeably Moreover,thefunctionofabiochemicaloscillatorneednotbe referred to as the oscillator. limited only to facilitate synchronization but can be extended to perform other functions such as providing timing infor- II. CALCIUMOSCILLATIONS mation for coordination among communication modules in Ca2+ oscillations are the fluctuations in [Ca2+]1 that arise a nanomachine or for scheduling channel access [6]. Hence, either from the entry of external calcium into a cell or from it would be a fair requirement of the biochemical oscillators the periodic release of intracellular calcium to the cytosol to operate and function as a separate unit within a nanoma- from internal stores [15]. The former is commonly found in chine but with the ability to communicate or pass timing excitablecellswhilethelatterinnon-excitablecells[16].Over information to other communication modules. The primary the years, several models have been developed to describe focus of our study is to address the former requirement, the mechanisms and dynamics of the Ca2+ oscillations. while the latter requirement is left as an open research For example, while some models consider Ca2+ oscillations challenge. with a constant Inositol 1,4,5-triphosphate (IP ) concentra- Based on the aforementioned issues and requirements, 3 tion [17]–[20], other modelsconsiderCa2+ oscillationswith a we propose a biochemical oscillator model by taking inspi- varying [IP ] [21]–[24]. ration from the widely investigated and studied biological 3 oscillatory phenomena called calcium (Ca2+) oscillations. The major mechanism and signaling principles of Ca2+ oscillationsfoundinastrocytes[23]isreviewedhereinabrief A bottom up approach is applied to model the biochemical manner. As shown in Fig. 1, in response to an external or oscillator as a separate functional unit within a nanomachine. agoniststimulation, the plasma membranetriggers the release Itemploysthreetypesofexcitatorymoleculesthatworkintan- of IP , which is a secondary messenger2, into the cytosol demtogenerateoscillatorypatternsintheconcentrationofthe 3 (denoted by ❶ ). The IP molecules propagate via diffusion molecule of interest — the molecule whose periodic increase 3 and bind to the bi-phasic IP receptors(IP Rs) located on the and decrease in concentration is considered to constitute the 3 3 membranes of Endoplasmic Reticulum (ER) (denoted by ❷ ). oscillator’s oscillations. The main contributions of this study The ER serves as an internal store of Ca2+ molecules. When are the design of a high frequency biochemical oscillator, the the IP Rs are bounded to by the IP molecules, they open investigationsto identify the parameters that affect the period 3 3 of the oscillations, and the identification of a key design unit. 1Squarebrackets denote theconcentration ofthechemical ormolecule. We state that a high frequency oscillator will allow faster 2Secondary messengers are molecules that relay information from a cell’s synchronization as the interval between the synchronization surfacetoothermolecules ororganelles thatresidewithinthecell. SHITIRIANDCHO:BIOCHEMICALOSCILLATORUSINGEXCITATORYMOLECULESFORNANONETWORKS 767 Fig.2. Ca2+ oscillations. up and allow the stored Ca2+ molecules to diffuse into the Fig.3. Internal signalingpathways oftheproposedoscillator. cytosol (denoted by ❸ ). This results in an increase in the cytosolic [Ca2+]. The increase in cytosolic [Ca2+] causes two actions.First,athigherconcentrationlevelsitbeginstoinhibit First, the oscillator is assumed to be equipped with an the IP Rs, thereby stopping the flow of Ca2+ from the ER to 3 internal store that holds B molecules; we call it the reserve the cytosol (denotedby ❹ ). Second, it stimulates the regener- unit(RU).ThesurfaceoftheRUisassumedtobecoveredwith ation of IP3 via PLCδ, also a secondary messenger (denoted receptorsthat are sensitive to both A and B molecules, called by ❺ ). Meanwhile, special structures called Ca2+ ATPase and the type AB sensitive receptors (tABRs), and with special Na+-Ca2+ exchanger regulate the cytosolic [Ca2+] to avoid suction-like receptors that are sensitive and activated only cell death caused by high levels of cytosolic [Ca2+]. The by B molecules when B’s concentration increases, called the Ca2+ATPase structures are found on the ER’s surface and the reserve unit pumps (RUPs). The aforementioned assumptions Na+-Ca2+ exchangerstructuresare foundon the cell’s surface are based on the internal Ca2+ stores and their receptors or membrane. Both, Ca2+ATPase and Na+-Ca2+ exchanger, (and pumps) [15] [16]. Then, we assumed the oscillator to be structures regulate the cytosolic [Ca2+] by pumping the Ca2+ a well-mixed3 closed unit [25], that is, molecules can neither intotheERandoutofthecell,respectivelyandarecommonly enter nor exit the oscillator unit. This assumption ensures referred to as pumps attributed to the pump-like action they two things. First, it ensures that unwanted noises, such as perform. A majority of the cytosolic Ca2+ is pumped into the the residual or lingering molecules from other processes that ER (denoted by ❻ ) and only a fraction of it is pumped out occurwithinthenanomachines,donotentertheoscillatorunit of the cell (not shown in the figure). The second action (❺ is and affect the oscillations. Second, it ensures the oscillator to responsible for increasing the cytosolic [Ca2+], while the first be a separate functional unit that can be embedded within action(❹ )combinedwiththepumpingactionsofCa2+(❻ )are a nanomachine. Next, the excitatory molecules are assumed responsible for reducing the cytosolic [Ca2+]. This periodic to have a lifetime and undergo biodegradation once they increaseand decreasein cytosolic[Ca2+] constitutein whatis accomplishtheirtask,therebyeliminatingthemselvesfromthe known as the Ca2+ oscillations shown in Fig. 2. oscillatorandminimizinganychancesofintra-oscillatorinter- Atthecellularlevel,theoscillationsof[Ca2+]ortheculmi- ferences. However, we consider the lifetime of B molecules nationof the oscillationsinto a wave is a formof information to be long enough so that the B oscillations can be sustained that can triggera variety of biologicalactivities dependingon overthelifetimeofthenanonetworkortheapplication.Finally, the type of cell in which they are produced. However, from while we show that a chemical process can regenerate A and a communication system point of view, the oscillations could theRU canstoreB, wesimplyassumeCtobesynthesizedby very well be interpreted as a clocking mechanism. As stated a chemical process [17]. For clarity and ease of readability, earlier, nanonetworks and its nanomachines are envisioned to we denote the concentration of B molecules in the RU and be largely based on molecules for their functionality,and this intra-oscillator space4 by [B] and [B] , respectively. RU INTRA indicates the importance and relevance of an oscillator that Fig. 3 depicts the proposed oscillator model along with employs molecules as well. Therefore, in the next section we the internal signaling pathways. The oscillator’s surface is introduce the oscillator model, describing the working mech- assumed to be covered with special receptors [1] that can anism along with elaborate discussions on the mathematical sense a beaconsignal(Bcon). Bcon isassumedto bespecific framework. moleculethatmayactasatriggertoinitializetheoscillator(s) ormaycarrythe informationrelatedtothesynchronizationof III. MODELOF THEBIOCHEMICAL OSCILLATOR oscillators. We will not elaborate on Bcon as it not the focus A. Working Principle of this paper and leave it aside for future work. When Bcon is sensed, the oscillator triggers A molecules. A molecules Theproposedoscillatormodelemploysthreedifferenttypes of excitatory molecules: A, B, and C to generate oscillatory patternsintheconcentrationofthemoleculeofinterest,which 3Theconcentrationsofthebiomoleculesarehomogeneousthroughout.This assumptionallowsustoobtainthesamemeasurementof[B]atanypointin is B. Before we describe the oscillator in detail, we provide theoscillator unit. the assumptions that are considered in this study. 4Volumeoftheoscillator excluding thereserveunit’svolume. 768 IEEETRANSACTIONSONNANOBIOSCIENCE,VOL.15,NO.7,OCTOBER2016 propagate via diffusion and bind to the tABRs present on the For our model, we require the flux balance equations pertain- surface of the RU. ing to two molecules, that of B and A. Now for the tABRs to able to open and close, we consider The changes in [B] are mainly governed by the INTRA each tABR to consist of three identical gate-like structures outward B flux through the tABRs and the inward B flux called sites. These gate-like structures have been observed in through the RUPs as shown in (1). J is the flux of B tABR biologicalreceptorssuchastheIPRs[17],[18].Thethreesites from the RU to the intra-oscillator space and J is the RUP aresiteA ,siteB ,andsiteB .Thesubscripts“actv” flux of B pumped from the intra-oscillator space into the RU. actv actv deactv and “deactv” denote the activation and deactivation functions Similarly, the balance equation for [A] can be formulated as performedbythesites.SiteA andsiteB areresponsible the difference between the maximum rate of A’s production actv actv fortheactivationofthetABRs,whilesiteBdeactvisresponsible stimulatedby B via C, JB+C→A, andthe loss rate of A, Jloss, forthe deactivationof the tABRs. Onlya single moleculecan as shown in (2). occupy each site at a time [26], e.g., a single A molecule d[B ] can bind to the Aactv site at a given time. When A molecules INTRA = JtABR − JRUP (1) dt bindto theA sites andBmoleculesbindto the B sites, actv actv d[A] only then the tABRs are activated and are said to be in the = JB+C→A− Jloss (2) dt open state. However, when B molecules occupy the B deactv site (occurs only when B concentration increases), then the Equations(1) and (2) are an over simplification and do not tABRs are deactivated and are said to be in the closed state. provide adequate details of the physiological processes that When the tABRs are in the open state, B molecules begin are involved in obtaining the changes in [B] and [A]. INTRA to diffuse out due to the concentration gradient ([B] is Hence, to add details of the physiological processes, the RU assumed to be higherthan [B] [17]).This results in an definitions and frameworks, hereon, provide the necessary INTRA increase in [B] . As [B] increases, B molecules mathematical formulations of the physiological processes. INTRA INTRA begin to occupy the third site, B , and the tABRs are We begin by modelling the dynamics of tABRs and as stated deactv forcedtogointotheclosedstate,therebystoppinganyfurther earlier, the tABRs can exist either in the open or closed state. outflow of B to the intra-oscillator space. Additionally, the Thus, we obtain the fraction of tABRs that are in the open RUPs begin to absorb B molecules back into the RU, aiding state (F ) [27] as shown by O in the reduction of [B] . By now, the concentration of INTRA N A and B molecules are decreasing, but in order to sustain F = O (3) O N the oscillations, B molecules have to be released again from T the RU. However, to release B molecules, the tABRs have where, N and N are the number of tABRs that are in O T to be opened and this means A molecules have to be regen- the open state and the total number of tABRs, respectively. erated as they are partly responsible for opening the tABRs. In addition, we stated that an open tABR allows B molecules Therefore,achemicalprocessinvolvingBandCmoleculesis to diffuse into the intra-oscillator space aided by the concen- considered for the regeneration or synthesis of A molecules. tration gradient ([B] −[B] ). However, for a tABR RU INTRA This chemical process is assumed to be initialized around the tobeintheopenstate,A andB sitesshouldbeactivated actv actv time when [B]INTRA is approaching its peak such that the by A and B molecules, respectively. To model this behavior, total time taken for A molecules to regenerate and propagate we incorporate the ion-channel formulation as described by to the tABRs is proportional to the time at which [B]INTRA Li and Rinzel [18]. Then the probability that Aactv and Bactv decreasesto a lower threshold.With the bindingof the regen- sites are in the activated state, such that an arbitrary tABR is erated A molecules (along with B the molecules of the intra- in the open state, is given by oscillator space) to the tABRs, B molecules are once again (cid:2) (cid:3)(cid:2) (cid:3) rTehleeaesendtirientporotcheessienstrac-uolmsciinllaatteorinstpoatcheeipnecrrieoadsiicngin[cBre]aIsNeTaRnAd. PB = [A][+A]K1 [B][IBN]TINRATR+AK2 (4) decreaseof[B] andarewhatconstitutestheoscillations INTRA where K and K are the dissociation constants5 of A and B of the proposed oscillator. 1 2 molecules from their respective sites. The first term signifies the activation of A sites and the second term signifies actv B. Mathematical Formulation the activation of the Bactv sites. Then, the product of (3) and (4) gives the overall probability of the tABRs that are In this section, we present the mathematical formulation in the open state. From these formulations, we can calculate of the oscillator model. For anydynamicsystems, such as the the flux of B from the RU, J , by taking into account oscillatorproposedinthispaper,thechangesinthesystemare tABR the concentration gradient ([B] −[B] ) that allows modelledusingordinarydifferentialequations(ODEs)andare RU INTRA B molecules to diffuse out with a rate r , the fraction of generally termed as the general flux balance equations. Flux, 1 tABRs that are in the open state, F , and the probability here,refersto the rate atwhich the moleculesflow,where the O of A and B sites that are in the activated state, P , motion of the flow is assisted by diffusion. In simple terms, actv actv B we need to identify the process involved in the generation 5Dissociationconstantdescribestheaffinitybetweenaprotein(e.g.,tABRs) and degradation of a molecule, and calculate the change in and a ligand (e.g., A). A small dissociation constant indicates a stronger the concentration of the molecule for every time instant. affinity. SHITIRIANDCHO:BIOCHEMICALOSCILLATORUSINGEXCITATORYMOLECULESFORNANONETWORKS 769 as shown by (5). The power, which is three, represents the Since A is dependent on B for its regeneration, we can number of sites per tABR. incorporate a factor α [17] to denote by how much the production of A depends on B. The value of α can range J =r (F )3(P )3([B] −[B] ) (5) tABR 1 O B RU INTRA between 0 and 1, with 0 denoting no dependenceat all and 1 denoting full dependence. Hence, we can rewrite (7) as Next, to quantify the flux of B molecules into the RU, aJsRU2P,[2w5e, ups.e2t9h3e]H. iHllilelqueaqtuioatniownitdhetshceriHbeilslcthoeeffibceiheanvtiovraluoef JB+C→A =r3[[BB]]INTRA++αKK4 (10) INTRA 4 molecules that have more than one binding site and a Hill’s Substituting(9),(6),(7),and(10)in(1)and(2),weobtainthe coefficientvaluegreaterthan1indicatespositivecooperativity complete flux balance equation for the oscillator as follows, among these binding sites. Positive cooperativity refers to the increase in the affinity of a biological substrate to other d[B]INTRA =η(r (F )3(P )3+r )([B] −[B] ) 1 O B 5 RU INTRA moleculesonceamoleculebindstoit.Equation(6)showsthe dt formulation of J using the Hill’s equation, [B]2 RUP − r INTRA (11) [B]2 2K32+[B]2INTRA JRUP =r2K2+[IBN]T2RA (6) d[A] =r [B]INTRA+αK4 −r [A] (12) 3 INTRA dt 3 [B] +K 4 INTRA 4 where,r isthemaximumrateatwhichBisabsorbedintothe 2 Now, given that the oscillator is assumed to be a closed RU and K is the dissociation constant of B molecules from 3 unit, we can apply the law of conservation of mass to the the RUPs. totalconcentrationof Bin theoscillator ([B] ). Therefore, TOT Thetermsonthe righthandside of(2)are merelytherates [B] is the sum of [B] and [B] , which implies TOT INTRA RU at which A is regenerated and degraded. As the regeneration that [B] can be calculated by RU of A is achieved by a chemical process involving B and C [B] −[B] molecules, this process can be modeled by [B] = TOT INTRA (13) RU η [B] JB+C→A =r3[B] INTR+AK (7) One final consideration to be made is to formulate the INTRA 4 rates at which the tABRs switch between the open states where, r3 and K4 are the maximum rate at which A can and the closed states. Such formulationis necessary to model be regenerated and the dissociation constant of B from C, the tABRs behavior as close as possible to the behavior respectively. Equation (7) is the Michaelis-Menten equation of the actual biological receptors and also to obtain N . O and is a special case of the Hill equation where the molecule Therefore,wemodelthestatetransitionsasatwostateMarkov has only one binding site. Meanwhile the degradationof A is process [27] [29, p. 286] proportional to a loss rate and is given by [17] J =r [A] (8) loss 4 where, r is the rate at which A is degraded. 4 Now,toprovideformaldetailstotheformulationsdiscussed where, Closed and Open refer to the closed state of tABRs so far (Equations (3) – (8)), first, we define and assign the when the B sites are bounded by B molecules and the deactv following biochemical molecules: type A, IP3; type B, Ca2+; open state of tABRs when B molecules are not bounded to andtypeC, PLCδ basedontheirresemblanceinfunctionality. them, respectively. The state transition parameters β and β o c Additionally, tABRs are defined by IP3Rs. Based on these are the rates at which the tABRs switch from the open states definitions, certain considerations have to be addressed in (5) to the closed states and vice versa, respectively. To determine and (8), and are discussed in the following paragraphs. the parameterswe apply the formulationdevelopedfor IP Rs 3 In(5),thedifferencesintheunitsofthefluxratesofBfrom by Jung and Shuai [26]. the RU and into the RU have not been taken into account. (cid:2) (cid:3) [A]+K The fact is the flux rate of B from the RU is measured in β =r K 1 (15) micromoles (μM) per liter RU per second, while the flux o 6 5 [A]+K6 rate of B into the RU is measured in micromoles per liter βc =r6[B]INTRA (16) intra-oscillatorper second.Thisis becausethe volumesof the where,r istherateatwhichBbindstothetABRs, K andK intra-oscillator space and the RU are not the same; volume 6 5 6 are the dissociation constants of B (inactivation) and A, of RU being smaller. To compensate for the differences in respectively. the volumes, a scaling factor η is introduced and is defined as the ratio of the intra-oscillator’s volume to the RU’s C. Operator, (cid:6) k volume [28, p. 23]. Moreover, as the RU is not considered To realize the characteristics befitting an oscillator, specifi- to be a closed unit, there will be a tiny amount of B leakage callytoattain highfrequencyofoscillations,we introducethe from the RU at a rate r . Therefore, (5) can be rewritten as, 5 concept of an operator ((cid:6) ). The operator is defined as the k J =η(r (F )3(P )3+r )([B] −[B] ) (9) factorbywhichavalueofcertainparameter,k,oftheoscillator tABR 1 O B 5 RU INTRA 770 IEEETRANSACTIONSONNANOBIOSCIENCE,VOL.15,NO.7,OCTOBER2016 can be varied to satisfy a criterion: attain high frequency of TABLEI oscillations while maintaining stable period of oscillations. PARAMETERVALUES We use the term “ variant ” to denote the parameter whose value is changed to the amount of modification made to it by the operator, all the while satisfying the criterion. Let Zk be the value of the kth parameter, then the variant value of the kth parameter, Zk , can be obtained by Var Zk = Zk(cid:6) , k ∈{(cid:7)} (17) Var k where, (cid:7) is the set of parameters that satisfy the criterion on being modified by the operator. Using Equation 17, we areabletoidentifytheparametersthatwillplayavitalrolein designingahighfrequencyoscillatorandthediscussionstothe same are providedin Section IV. To highlightthe importance of faster oscillations for faster synchronization, let us outline two case examples. In both cases, we consider synchronization using a pulse- coupling like scheme [30], where the oscillators sense the oscillation peaks of the other nodes through the synchro- nization signals such as Bcon, adjust their peaks accord- ingly, and over a period of time converge towards peaking at the same times. In the first case, a reference oscillator sends out Bcon signals at the beginning of every period, a period being made up of the time between two consecutive peaks. Now, if the reference oscillator has a period of, say, 2 seconds, it means that the Bcon signal can be sent after every2seconds.Butiftheoscillatorhasalargerperiodof,say, 5 seconds, then the Bcon signal can be sent out after every 5 seconds. Intuitively, notwithstanding the stochastic nature of the molecular channel, the Bcon sent out from the faster oscillator (period = 2 seconds) will arrive at a specific destination much earlier or will have a higher probability of reaching earlier than that sent out from the slower oscillator (period= 5 seconds).In the secondcase, pulse-couplingwith Fig. 4. Illustration of two consecutive periods having certain a two-way exchange mechanism, where a reference oscillator periods Xi and Xi+1. sends out the Bcon signal only when it receives a response signal from the other oscillators is considered. Such mech- anism avoids redundant signaling as the reference oscillator thevaluesoftheparametersusedinthesimulations.Aparame- needs to send out the Bcon signal only if the response signal ter’svalueissaidtobetheoptimumvalueifthatvalueresults arrives before the beginning of an immediate period; if not, inthesmallestmeanandvarianceintheperiodofoscillations. Bconisnottransmitted.Again,intuitively,iftheoscillatorhas Period of oscillation (Xi,i =1,2,3...n−1) refers to the a smaller period (e.g., 2 seconds), then the Bcon signal can time between consecutive peaks, where “n” is the number of be sent out at the immediate peak, which will be at a shorter peaks observedfor t seconds. If X is a random variable, then time, allowing Bcon to reach at a destination much earlier. the varianceof the peri(cid:4)odcan be obtainedas σ =(Xi − X¯)2, Besidesfacilitatingfastersynchronization,ahighfrequency withmeanperiod X¯ = X /n−1asshowninFig.4.Unless i oscillator has another benefit and that is it will allow higher specified or denoted by ∗, the values of the parameters listed data rates. For instance, when a sender has data to send, it are the values used in [17], which we refer to as the standard has to encode the message through, say, the concentration values (Zk ). Std or frequency of the message molecule(s). The rate, at which We define (cid:7) as r , r , r ,K , and α, thus, k ∈ {r , r , r , 1 2 3 3 1 2 3 the data are encoded, regardless of the encoding scheme, is K ,andα}and Zk isreplacedby Zk .Theaforementionedk 3 Std proportionalto the oscillator’sfrequency— a faster oscillator parameterswere identifiedthroughextensivesimulationssuch meanshigherdatarate.Andtherefore,havingahighfrequency that the criterion is satisfied. In addition, the following bases oscillator can benefit a nanonetwork in two significant ways. werealso usedindeterminingthe parameters.Theparameters r , r , and K were determined so that the resultant B oscil- 1 2 3 IV. SIMULATIONRESULTSAND lations are bounded within reasonable [A]; while r and α 3 ANALYSIS OF PARAMETERS were determined to maintain the balance in (9) and (10). The simulations were carried out on a time scale Thus, from these bases, extensive simulations were carried of 100 seconds with a time step of 0.01 seconds. Table I lists out to determine the optimum values of the k parameters SHITIRIANDCHO:BIOCHEMICALOSCILLATORUSINGEXCITATORYMOLECULESFORNANONETWORKS 771 TABLEII PARAMETERSWITHTHEIRSTANDARD,OPERATOR, ANDVARIANTVALUES theiroptimumvaluescauseanincrease,whichisnotdesirable, in the mean and variance. Therefore, one must be careful in choosing these values. In Fig. 5(c) and Fig. 5(d), the graphs followadecreasingslopeforthemeanandanincreasingslope for the variances. For both of the figures, a trade-offbehavior isobserved,thatis, forlowervaluesofmean,largervariances are obtained and vice versa, suggesting that the optimum values of the mean and variance for r and k lie in and 3 3 around the point of intersection. In Fig. 5(e), increase in α cause the graphs to decrease. Note that the graphs after 0.6 is actually a slope but with very small differences (≈4 μs), and it gives the impressions of a straight line. Therefore, a minimum of ≈60% dependencerate of A must be allowed to achieve negligible variance and uniform periods. It should be notedthatinFig.5(a)–(e),therangesonthex-axesarelimited only to those values beyond which the oscillations die out or the mean and variance are too large to consider the oscillator suitable for a communication nanonetwork. Table II lists the operator value associated to each standard value that results in the corresponding variant values. To show how much of a difference was obtained in the frequencies of oscillations using the variant values of the proposed oscillator and the standard values of the standard oscillator [17], we plot the oscillations for each oscillator in Fig. 6. As expected, in Fig. 6(a), a higher frequency of oscillationwithameanperiodof2.0148secondswasobtained for the variant values, while in Fig. 6(b) a lower frequency of oscillations with a mean period of 17.71 seconds was obtained for the standard values. The oscillations in Fig. 6(a) were achieved by using the concept of the operator and the associated criterion. One may notice the difference in the “optimum” values of the variant parameters used (Table I) to those obtained in Fig. 4 (a)–(e). That is because we numerically evaluated and chose the variant values that were more suitable in gen- erating the resultant oscillations shown in Fig. 6(a). However, the variant values used satisfy the criterion and reside well within the range of the optimum values, and therefore, has Fig. 5. Mean and variance of B oscillations for different values of (a) r1 insignificantornoeffectontheoutcome.Infact,asanexample (b)r2 (c)r3 (d) K3 and(e)α. to justify that difference, we plot the mean and variance as shown(cid:4)in Fig. 5 (a)–(e). Hereon, unless(cid:4)specified, mean versus two parameters, r1 and r2, simultaneously. r1 and r2 refers to X¯ /N and variance refers to σ /N; where were specifically chosen due to their association to the RU. j j j =1,2,3...N and N is the number of simulation runs. Atrade-offcanbeobservedbetweenthemeanandthevariance In Fig. 5(a) and Fig. 5(b) the graphs of the mean and (Fig. 7 (a), (b)) but the key observation to note here is that variance follow a valley-like shape, suggesting the optimum the values of r and r (ref. Table 1) results in a very low 1 2 values of r and r lie in and around the valleys. Both these variance and a relatively decent mean (denoted by the white 1 2 parameters are associated with the RU and deviations from cross). Another observation to note is that a large number of 772 IEEETRANSACTIONSONNANOBIOSCIENCE,VOL.15,NO.7,OCTOBER2016 Fig.8. Meanandvariance fordifferent initial concentrations ofA. Fig.9. Variancefordifferent numberofreceptors, NTOT. Fig.6. Boscillations with(a)variant values (b)standardvalues. whichmolecularcommunicationoccurs,itisappropriatetobe termedashigh.TheinitialconcentrationofAfortheproposed oscillatorissetto2μMasthisvalueaccountsforthesmallest mean as shown in Fig. 8, while that of the standard oscillator is set to 0.5 μM [17]. On careful examination of the values listed in Table II, the operator’s value associated to parameters r and r are 1 2 quite large, indicating that the tABRs and RUPs should be designed to release and absorb B molecules, respectively, in larger quantities or be able to perform their functions at a faster rate. Either way their design will play a key role in satisfying the criterion mentioned in Section III. Another parameterthat couldplaya vitalrole in maintainingthe crite- rion is the number of tABRs, N . A higher N means TOT TOT that more tABRs will be in the open state, allowing more B molecules to be released into the intra-oscillator space at a given instant of time. This, in a way, also aides in generating smoothoscillationcurves,whichotherwisewouldbejaggedor spiky caused by the random nature in which the B molecules are released. Therefore, with higher N , [B] will TOT INTRA increase faster and it results in a proportionate change in the regeneration of A as well as the inhibition of tABRs. As a result, shorter periods can be obtained. Furthermore, a high N is desirable to ensure that the variances are very small TOT Fig. 7. (a) Mean and (b) variance obtained by varying the ornegligible.AsshowninFig.9,thevariancetendstobecome parameters r1 andr2 simultaneously. smaller as NTOT is increased incrementally from 1×105 to 10×105 and a gentle slope can be observed when N is TOT equal to or greater than 5×105. the values of variance reside in the lows, and therefore most From these observations, we can infer that the parameters of the combinations of r and r will lead to a minimum in N , r , and r have great influences on the oscillator’s 1 2 TOT 1 2 the variance. This may allow flexibility in choosing any of behaviorandinterestinglyallthreeareassociatedwiththeRU. the combinations of r and r within the lows as long as the This suggests the design of RU to be a key role in realizing a 1 2 criterion is not violated. stable high frequency biochemical oscillator. That being said, In terms of magnitude, the high frequency of oscillations it remains an open research challenge to identify methods on for biochemical oscillators will be nowhere close to those how such reserve unit(s) can be physically designed or to of traditional oscillators. However, considering the speed at be more specific how the operator’s value can be physically SHITIRIANDCHO:BIOCHEMICALOSCILLATORUSINGEXCITATORYMOLECULESFORNANONETWORKS 773 realized.Notwithstandingthechallenge,westronglybelievein [2] S. Hiyama, T. Nakano, A. W. Eckford, and T. Haraguchi, “Molecular the viability of the realization of such structuresin the future, communication,”inProc.NSTINanotechnol.Conf.TradeShow,vol.3. 2005,pp.391–394. ifnotnow,aswehavetakenutmostcareinusingassumptions [3] R.A.Freitas,Jr.,Nanomedicine:BasicCapabilities,vol.1.BocaRaton, thatarerealisticandgenerallyaccepted.Finally,theinfluences FL,USA:CRCPress,2012,ch.2. of the remaining parametersr , K , and α on the oscillator’s [4] T. Suda, M. Moore, T. Nakano, R. Egashira, and A. Enomoto, 3 3 “Exploratory research on molecular communication between nanoma- behavior, though minimal, cannot be ignored as they provide chines,”Nat.Comput.,vol.25,pp.1–30,Jan.2005. stability to the system when r1 and r2 are varied. [5] I. F. Akyildiz, F. Fekri, R. Sivakumar, C. Forest, and B. K. Hammer, “Monaco: Fundamentals of molecular nano-communication networks,” V. CONCLUSION IEEEWireless Commun.,vol.19,no.5,pp.12–18,May2012. [6] T. Nakano, T. Suda, Y. Okaie, M. J. Moore, and A. V. Vasilakos, In this paper, a biochemical oscillator for nanonetworks “Molecular communication amongbiological nanomachines: Alayered was proposed. The major difference of the study from other architecture and research issues,” IEEE Trans. Nanobiosci., vol. 13, literature was the use of excitatory molecules to generate no.3,pp.169–197,Sep.2014. [7] M. J. Moore and T. Nakano, “Oscillation and synchronization of the oscillations. Besides this, we also designed the oscillator molecular machines by the diffusion of inhibitory molecules,” IEEE as a separate functional unit that can be embedded within Trans.Nanotechnol., vol.12,no.4,pp.601–608, Jul.2013. a nanomachine. A detailed mathematical formulation of the [8] U.AkgülandB.Canberk,“Aninterference-freeandsimultaneousmole- oscillatorwasprovided,wherein,theconceptofan“operator” culartransmissionmodelformulti-usernanonetworks,”NanoCommun. Netw.,vol.5,no.4,pp.83–96,2014. was introduced to show that a high frequency biochemical [9] U. Alon, An Introduction to Systems Biology: Design Principles of oscillator can be obtained. Biological Circuits. BocaRaton, FL,USA:CRCPress,2006. ExtensivesimulationswerecarriedoutinMATLABtostudy [10] M. B. Miller and B. L. Bassler, “Quorum sensing in bacteria,” Annu. Rev.Microbiol., vol.55,pp.165–199,Oct.2001. the behavior of the oscillator when certain parameters were [11] A.Goldbeter, “Computational approaches tocellular rhythms,”Nature, changed.The quantitativeanalysis of the oscillator’s behavior vol.420,pp.238–245,Nov.2002. led to the identification of five parameters, viz r , r , r , [12] M.B.ElowitzandS.Leibler, “Asynthetic oscillatory networkoftran- 1 2 3 K , and α, whose modification by the operator (towards an scriptional regulators,” Lett.Nature,vol.403,pp.335–338,Nov.1999. 3 [13] J.Garcia-Ojalvo, M.B.Elowitz, andS.H.Strogatz, “Modeling asyn- optimum value) resulted in obtaining high frequency stable thetic multicellular clock: Repressilators coupled by quorum sensing,” oscillations. In addition to these parameters, N was also Proc.Nat.Acad.Sci.USA,vol.101,no.30,pp.10955–10960, 2004. TOT [14] M. Toiya, H. O. Gonzalez-Ochoa, V. K. Vanag, S. Fraden, and identifiedtoplayavitalroleinobtaininghighfrequencystable I.R.Epstein,“Synchronization ofchemical micro-oscillators,” J.Phys. oscillations.Thestudyalsoincludesqualitativeanalysisofthe Chem.,vol.1,no.8,pp.1241–1246, 2010. parametersand identifiedthe parameters N ,r , andr , to [15] M.J.Berridge,“Calciumoscillations,”J.Biol.Chem.,vol.265,no.17, TOT 1 2 p.9583,1990. have great influence on the oscillator’s behavior. Interestingly [16] G. Dupont, L. Combettes, G. S. Bird, and J. Putney, “Calcium oscil- allthreeparametersareassociatedwiththeRU.Thisindicates lations,” Cold Spring Harbor Perspect. Biol., vol. 3, p. a004226, and suggests the importance of the RU’s design in realizing Mar.2011. [17] G. W. De Young and J. Keizer, “A single-pool inositol the high frequency stable biochemical oscillator. 1,4,5-trisphosphate-receptor-based model for agonist-stimulated The proposed oscillator is envisioned not only to facilitate oscillations in Ca2+ concentration,” Nature, vol. 89, pp. 9895–9899, synchronization but also to facilitate other functions such as Oct.1992. providing timing information for scheduling channels access [18] Y.X.LiandJ.Rinzel,“EquationsforInsP3 receptor-mediated[Ca2+]i oscillations derived from adetailed kinetic model: A Hodgkin–Huxley and decoding the signals or for coordinating other commu- likeformalism,”J.Theor.Biol.,vol.166,no.4,pp.461–473, 1994. nication modules in a nanomachine. For these reasons, and [19] M.Falcke, “Reading thepatterns inlivingcells—ThephysicsofCa2+ besides making the oscillator free from unwanted external signaling,” Adv.Phys.,vol.53,no.3,pp.255–440, 2004. [20] Z. Shuai, L. Bing, Z. Shaoqun, and C. Shangbin, “Simulation of noises, the oscillator is envisioned as a separate functional spontaneous Ca2+ oscillations in astrocytes mediated by voltage- unit that can be embeddedwithin a nanomachine.Apart from gated calcium channels,” Biophys. J., vol. 97, no. 9, pp. 2429–2437, theaforementionedfunctions,ananonetworkcanbenefitfrom 2009. [21] T. Meyer and L. Stryer, “Molecular model for receptor-stimulated a high frequency oscillator in terms of faster data rates. calcium spiking,” Proc. Nat. Acad. Sci. USA, vol. 85, pp. 5051–5055, In conclusion, the rationale for undertaking this study is Jul.1988. to provide a physical clocking mechanism for nanomachines [22] U. Kummer, B. Krajnc, J. Pahle, A. K. Green, C. J. Dixon, and M. Marhl, “Transition from stochastic to deterministic behavior in that will, in general, facilitate synchronization among the calcium oscillations,” Biophys.J.,vol.89,no.3,pp.1603–1611,2005. nanomachines in a nanonetwork. Being minuscule, nanoma- [23] T.Höfer,L.Venance,andC.Giaume,“Controlandplasticityofintercel- chinesarelimitedinfunctionalityandpower.Thus,toperform lularcalcium waves inastrocytes: Amodelingapproach,” J.Neurosci., vol.22,no.12,pp.4850–4859, 2002. complex and large operations would require them to coordi- [24] M. De Pittà, M. Goldberg, V. Volman, H. Berry, and E. Ben-Jacob, nate or synchronize with other nanomachines, for example, “Glutamate regulation of calcium and IP3 oscillating and pulsating releasing drugs to a targeted site. This makes synchrony dynamics in astrocytes,” J. Biol. Phys., vol. 35, no. 4, pp. 383–411, 2009. even more crucial for nanonetworks and remains an open [25] J.KeenerandJ.Sneyd,Mathematical PhysiologyI,2nded.NewYork, research challenge. Therefore, in the future we will study the NY,USA:Springer, 2008. design of synchronization schemes for nanonetworks, whose [26] J.-W. Shuai and P. Jung, “Stochastic properties of Ca2+ release of nanomachines are equipped with the proposed oscillator. inositol1,4,5-trisphosphatereceptorclusters,”Biophys.J.,vol.83,no.1, pp.87–97,2002. REFERENCES [27] J. W. Shuai and P. Jung, “Optimal ion channel clustering for intracel- lular calcium signaling,” Proc. Nat. Acad. Sci. USA, vol. 100, no. 2, [1] I. F. Akyildiz, F. Brunetti, and C. Blázquez, “Nanonetworks: pp.506–510, 2003. A new communication paradigm,” Comput. Netw., vol. 52, no. 12, [28] J. Sneyd, Tutorials in Mathematical Biosciences II Mathematical. pp.2260–2279, Aug.2008. NewYork,NY,USA:Springer, 2005. 774 IEEETRANSACTIONSONNANOBIOSCIENCE,VOL.15,NO.7,OCTOBER2016 [29] C.P.Fall,Comput.CellBiol.,vol.20,2004. Ho-ShinCho(S’92–M’99)receivedtheB.S.,M.S., [30] R. E. Mirollo and S. H. Strogatz, “Synchronization of pulse- andPh.D.degreesinelectricalengineeringfromthe coupled biological oscillators,” SIAM J. Appl. Math., vol. 50, no. 6, KoreanAdvancedInstitute ofScienceandTechnol- pp.1645–1662, 1990. ogy(KAIST)in1992,1994,and1999,respectively. From 1999 to 2000, he was a Senior Member of theResearchStaffwiththeElectronicsandTelecom- munications Research Institute (ETRI), where he EthungshanShitiri(S’16)receivedtheB.E.degree was involved in developing a base station system inelectronics andcommunication engineering from forIMT-2000.From2001to2002,hewasafaculty Thanthai Periyar Government Institute of Tech- member at the School of Electronics, Telecommu- nology, India, in 2010 and the M. Tech. degree nications, andComputer Engineering, Hankuk Avi- in communication systems from Christ University, ation University. In2003, hejoined the faculty ofthe School ofElectronics India, in 2013. He is currently working toward the Engineering, Kyungpook National University, South Korea, and now is a Ph.D.degreeatMobileCommunicationLaboratory, Professor. His current research interests include mobile communication— Kyungpook National University, South Korea. His IoT, 5G cellular; underwater communication—CDMA, SON, NOMA; and currentresearchinterestsincludemolecularcommu- molecular communication—synchronization andmultiplexing. nication,synchronizationandmediumaccesscontrol Prof. Cho was awarded the rising researcher fellowship of NRF in 1998, forwireless networks,suchasnanonetworks. andisamemberofIEEE,IEICE,IEEK,KICS,andASK.