Research article Special Issues

Fully Bayesian analysis of allele-specific RNA-seq data

  • Diploid organisms have two copies of each gene, called alleles, that can be separately transcribed. The RNA abundance associated to any particular allele is known as allele-specific expression (ASE). When two alleles have polymorphisms in transcribed regions, ASE can be studied using RNA-seq read count data. ASE has characteristics different from the regular RNA-seq expression: ASE cannot be assessed for every gene, measures of ASE can be biased towards one of the alleles (reference allele), and ASE provides two measures of expression for a single gene for each biological samples with leads to additional complications for single-gene models. We present statistical methods for modeling ASE and detecting genes with differential allelic expression. We propose a hierarchical, overdispersed, count regression model to deal with ASE counts. The model accommodates gene-specific overdispersion, has an internal measure of the reference allele bias, and uses random effects to model the gene-specific regression parameters. Fully Bayesian inference is obtained using the fbseq package that implements a parallel strategy to make the computational times reasonable. Simulation and real data analysis suggest the proposed model is a practical and powerful tool for the study of differential ASE.

    Citation: Ignacio Alvarez-Castro, Jarad Niemi. Fully Bayesian analysis of allele-specific RNA-seq data[J]. Mathematical Biosciences and Engineering, 2019, 16(6): 7751-7770. doi: 10.3934/mbe.2019389

    Related Papers:

    [1] Eduardo Ibarguen-Mondragon, Lourdes Esteva, Leslie Chávez-Galán . A mathematical model for cellular immunology of tuberculosis. Mathematical Biosciences and Engineering, 2011, 8(4): 973-986. doi: 10.3934/mbe.2011.8.973
    [2] Maysaa Al Qurashi, Saima Rashid, Fahd Jarad . A computational study of a stochastic fractal-fractional hepatitis B virus infection incorporating delayed immune reactions via the exponential decay. Mathematical Biosciences and Engineering, 2022, 19(12): 12950-12980. doi: 10.3934/mbe.2022605
    [3] Bruno Buonomo, Marianna Cerasuolo . The effect of time delay in plant--pathogen interactions with host demography. Mathematical Biosciences and Engineering, 2015, 12(3): 473-490. doi: 10.3934/mbe.2015.12.473
    [4] Maria Vittoria Barbarossa, Christina Kuttler, Jonathan Zinsl . Delay equations modeling the effects of phase-specific drugs and immunotherapy on proliferating tumor cells. Mathematical Biosciences and Engineering, 2012, 9(2): 241-257. doi: 10.3934/mbe.2012.9.241
    [5] Eduardo Ibargüen-Mondragón, Sandra P. Hidalgo-Bonilla, Miller Cerón Gómez . On macrophage response to primary Mycobacterium tuberculosis in humans. Mathematical Biosciences and Engineering, 2025, 22(9): 2506-2525. doi: 10.3934/mbe.2025092
    [6] Urszula Foryś, Jan Poleszczuk . A delay-differential equation model of HIV related cancer--immune system dynamics. Mathematical Biosciences and Engineering, 2011, 8(2): 627-641. doi: 10.3934/mbe.2011.8.627
    [7] Mohammed Meziane, Ali Moussaoui, Vitaly Volpert . On a two-strain epidemic model involving delay equations. Mathematical Biosciences and Engineering, 2023, 20(12): 20683-20711. doi: 10.3934/mbe.2023915
    [8] Jordi Ripoll, Jordi Font . Numerical approach to an age-structured Lotka-Volterra model. Mathematical Biosciences and Engineering, 2023, 20(9): 15603-15622. doi: 10.3934/mbe.2023696
    [9] Suqi Ma . Low viral persistence of an immunological model. Mathematical Biosciences and Engineering, 2012, 9(4): 809-817. doi: 10.3934/mbe.2012.9.809
    [10] Beryl Musundi . An immuno-epidemiological model linking between-host and within-host dynamics of cholera. Mathematical Biosciences and Engineering, 2023, 20(9): 16015-16032. doi: 10.3934/mbe.2023714
  • Diploid organisms have two copies of each gene, called alleles, that can be separately transcribed. The RNA abundance associated to any particular allele is known as allele-specific expression (ASE). When two alleles have polymorphisms in transcribed regions, ASE can be studied using RNA-seq read count data. ASE has characteristics different from the regular RNA-seq expression: ASE cannot be assessed for every gene, measures of ASE can be biased towards one of the alleles (reference allele), and ASE provides two measures of expression for a single gene for each biological samples with leads to additional complications for single-gene models. We present statistical methods for modeling ASE and detecting genes with differential allelic expression. We propose a hierarchical, overdispersed, count regression model to deal with ASE counts. The model accommodates gene-specific overdispersion, has an internal measure of the reference allele bias, and uses random effects to model the gene-specific regression parameters. Fully Bayesian inference is obtained using the fbseq package that implements a parallel strategy to make the computational times reasonable. Simulation and real data analysis suggest the proposed model is a practical and powerful tool for the study of differential ASE.


    1. Introduction

    The increasing penetration of renewable generation has been motivated by the necessity of reducing the current dependence on non-renewable energy sources. Solar and wind power systems are two of the most widespread forms of electricity generation based on renewable energy resources. The modular nature of these technologies makes them a suitable option for the deployment of small scale generation connected at distribution level, namely as distributed generation (DG).

    DG can be used for supporting voltage, reducing losses, providing backup power and ancillary services, or deferring distribution system upgrade [1,2,3]. However, the intermittent nature of some renewable resources (e.g. solar and wind) complicates their integration into the grid as they cannot be properly dispatched. In general, renewable generation is permitted to produce as much power as the resource availability allows; this approach limits the integration of renewable DG since unacceptable operating conditions may appear under large penetration scenarios (e.g. system voltages above accepted limits).

    Energy storage (ES) capabilities can be used to improve network operation and compensate the intermittency of renewable generation [4,5,6,7]. Through the aggregation of renewable DG and ES, it is possible to constitute an entity, known as a virtual power plant (VPP), which can operate as a single dispatchable power plant [8,9]. As part of a VPP, ES can be charged when there is a generation excess and provide support to meet the dispatched power when DG generation decreases.

    This paper presents the implementation of a custom-made VPP model for power flow calculations in OpenDSS [10]. The objective is to present a model based on the storage object available in OpenDSS, adequate for distribution system studies and quasi-static (time-driven) simulations. The VPP model has been compiled within the OpenDSS COM DLL and can be used in simulations driven from other software platforms, e.g. MATLAB.

    The paper is organized as follows. Section 2 details the configuration and operation of the implemented model. A case study aimed at demonstrating the model behavior is presented in Section 3. The optimization of the hourly day-ahead dispatch using a genetic algorithm (GA) is explored in Section 4. The main conclusions drawn from the paper, as well as future work, are summarized in Section 5.


    2. VPP Model for Power Flow Calculations


    2.1. VPP definition

    A VPP can be defined as a flexible representation of a portfolio of distributed energy resources (DERs), in which a single operating profile can be created from a composite of parameters characterizing each DER [11]. Through the aggregation of DERs it is possible for the VPP to operate as a system-connected power plant and have access to wholesale energy markets and provide support services for the transmission system management.

    VPPs can be classified as commercial VPP (CVPP) or technical VPP (TVPP) [11]. A number of DERs can participate in the energy market as a single CVPP, while the TVPP models the characteristics of a system that contains an aggregation of DERs, and acts as a single power plant. The CVPP can represent resources from different locations, whereas all DERs in a TVPP belong to the same geographical location.

    The main elements that compose a VPP are: (ⅰ) energy production units, (ⅱ) energy storage units, and (ⅲ) flexible loads [12]. The VPP must exert control over the different elements in order to coordinate and optimize their operation. VPP control is conducted via a communication system that connects all the DERs in the VPP [8,13].

    A VPP can also be defined as an aggregation of power production from a group of grid-connected DERs operated by a centralized controller [9]. According to this definition, the combination of DG and ES with a coordinated control can also be considered as a VPP, since it allows them to operate as a single dispatchable power plant. Although exporting energy to the transmission system is still possible (by reversing the power flow through the substation transformer), the VPP's main task will be to provide support to the distribution network as dispatchable DG. The model presented in this section corresponds to this definition of VPP and only considers the aggregation of energy storage and renewable generation, either photovoltaic (PV) or wind.


    2.2. VPP configuration and components

    The general configuration of the VPP model is shown in Figure 1; it is composed by a generation unit (GU), an ES unit, and an inverter. In this configuration the GU, which can be either solar or wind generation, and the ES are connected at DC level. The VPP will be referred to as Solar VPP or Wind VPP depending on the nature of the renewable resource used in the GU. The VPP is capable of injecting or absorbing reactive power through the inverter. The required parameters for each component are listed in Table 1.

    Figure 1. Basic configuration of the virtual power plant.
    Table 1. Parameters of a virtual power plant.
    PV Generator Energy Storage Unit
    •Rated power (kW) •Rated power (kW)
    •Irradiance (kW/m 2) •Rated energy (kWh)
    •Power vs. Temperature curve •Minimum stored energy (kWh)
    •Panel temperature (℃) •Efficiency curve
    •Minimum power (%) •Minimum power (%)
    •Power curtailment (%)
    Wind Generator Inverter
    •Rated power (kW) •Rated power (kVA)
    •Wind speed (m/s) •Rated voltage (kV)
    •Hub height (m) •Power factor
    •Cut-in speed (m/s) •Efficiency curve
    •Cut-out speed (m/s) •Minimum power (%)
    •Power vs. Wind speed curve •No-load power (%)
    •Wind Reference height (m) •Number of phases
    •Altitude over sea-level (m) •Connection (delta/wye)
    •Roughness length in wind direction (m) •ZIP coefficients
    •Ambient temperature (℃)
    •Power curtailment (%)
     | Show Table
    DownLoad: CSV

    PV Generator: The power produced by a PV generator is a function of the irradiance, panel temperature, and rated power at a selected temperature and an irradiance of 1.0 kW/m2[10]. It is assumed that the model can operate at maximum power or reduce the generated power according to the power curtailment parameter. Furthermore, the PV generator will only inject power if the generated power is above the predetermined minimum power, which is defined as a percentage of the rated power.

    Wind Generator: The output power produced by the wind generator is calculated following the procedure presented in [14]. The wind generator is defined by its rated power, hub height, cut-in and cut-out speed, and power vs. wind curve. The wind speed is defined for a reference height and the corresponding height and air density corrections are applied. As the PV generator, the wind generator also has a minimum generated power, which is defined as a percentage of the rated power.

    ES Unit: It is a generic technology-independent model defined by its rated power, rated energy, minimum stored energy, efficiency (which is a function of ES loading) and a minimum charge/discharge power (defined as a percentage of the ES rated power); see [10,15]. Moreover, the model automatically sets a maximum charge/discharge power based on the stored energy and the simulation time-step.

    Inverter: The inverter acts as a bridge between the DC network, to which the GU and ES unit are connected, and the AC distribution system; it allows a bidirectional power flow between the VPP and the grid. The inverter only absorbs/injects power if the power on the AC-side is greater than the predefined minimum power (defined as a percentage of the inverter rated power). The efficiency curve is defined as a function of the power output (in kVA) with respect to the inverter rated power. The absorbed power during the idle conditions (i.e. no active power output) is determined by the No-load power value. The idle power consumed by the VPP depends on the passive elements (e.g. filters and transformers) and the necessary power to operate the DC network; this power is defined in the model as a percentage of the inverter rated power and modeled as a ZIP load. The ZIP load model is constructed as a combination of constant power, constant current and constant impedance load models [16]; ZIP coefficients are defined in the following order:

    ZIPcoefficients=[CZ,CI,CP] (1)

    where CZ is the weight assigned to the constant impedance load model, CI is assigned to the constant current model, and CP is the weight that represents the portion of the load modeled as constant power.

    The inverter is also responsible for providing the required reactive power, being the maximum reactive power that the inverter can inject/absorb determined by the user-defined power factor (PF) parameter:

    Qmax=kVARATED1PF2 (2)

    The provided reactive power will have an effect on the inverter efficiency; this effect will be included, using the same approach presented in [17], as follows:

    η(s,pf)=f(s)α(pf) (3)
    s=kVAkVARATED (4)
    α(pf)=0.98+0.02pf (5)

    where Qmax is the maximum reactive power and kVARATED is the inverter rated power. η is the efficiency curve for any load level, s, and any output power factor, pf, f (s) is the efficiency curve for a unity power factor, which depends on the load level (measured at the AC-side), and α is a scale factor that depends on the power factor.


    2.3. VPP operation

    The VPP model has been developed for time-driven simulations and takes advantage of OpenDSS built-in capabilities (i.e. loadshapes, linear interpolation, complex number operations, etc.). Two operation modes have been defined for the VPP: (ⅰ) Follow mode; (ⅱ) Dispatch mode.


    2.3.1. Follow operation mode

    Under this operation mode the ES unit and the GU behave according to the curve shapes assigned to each component. The PV generator requires curves describing the solar irradiance and panel temperature for the evaluation period, while the wind generator follows the curves that define the wind speed (measured at a reference height) and the ambient temperature. On the other hand, ES unit operates according to a curve that defines the charge and discharge states, as well as the absorbed and injected power. The reactive power provided by the inverter is equal to the maximum reactive power times the corresponding multiplier stored in a user-defined curve shape.

    The total VPP power output is not controlled and depends on the power resulting from the interaction between the GU and the ES unit. The model follows the same convention as the OpenDSS Storage object: Negative values correspond to power absorption, while power injection is defined by positive values. The VPP powers are calculated according to the following equations (see Figure 1):

    PIN=PDG+PES (6)
    POUT={ηPIN,POUT0PINη,POUT<0 (7)

    where PDG is the power produced by the GU, PES is the power provided by the ES unit, PIN is the power measured at the inverter DC-side, POUT is the active power measured at the inverter AC-side, and η is the inverter efficiency, calculated according to Equations (3) through (5).


    2.3.2. Dispatch operation mode

    Under this operation mode, the VPP is required to follow a dispatch curve that defines the active and reactive power output for each time-step during the evaluation period. The objective is to fully utilize the capacity of the renewable resource, so ES operation will be subjected to the lack or excess of renewable generation. Additionally, the VPP can draw power from the grid in order to charge the ES unit. The complete operation strategy is depicted in Figure 2; note that Equations (6) and (7) remain valid.

    Figure 2. Dispatch operation mode.

    In Figure 2 EES is the current stored energy in the ES unit, EES_RATED is the ES rated energy, EES_MIN is the minimum stored energy, PES_MAX is the maximum power the ES unit can absorb/inject based on the current stored energy, and pc is the power curtailment parameter.

    The following aspects must be taken into account for the Dispatch operation mode:

    POUT is calculated as the inverter rated power times the corresponding multiplier stored in a user-defined curve shape.

    • The reactive power provided by the inverter (QOUT) is equal to the maximum reactive power times the corresponding multiplier stored in a user-defined curve shape.

    PIN is calculated from Equation (7).

    • The power produced by the GU is defined by the curves describing the renewable resource (irradiance or wind speed) and temperature (panel or ambient) for the evaluation period.

    • The ES unit main task is to compensate for the lack or excess of renewable generation. ES power will be adjusted so PIN matches the power required to meet the VPP dispatch.

    • By defining a negative multiplier, the VPP can be forced to absorb active power, which will be used to charge the ES unit with power coming from the distribution network.

    • VPP operation is meant to follow the power dispatch as accurately as possible; therefore, renewable generation may be curtailed in order to prevent any deviation from the dispatched power. For example, assume the ES unit is fully charged and the available power is greater than the required power to meet the VPP dispatch; under this condition, the GU production will be reduced in order to match the dispatched power. However, power curtailment is not a binary decision, that is, the amount of power to be curtailed is defined by the power curtailment parameter.

    • If the available power in the VPP (including ES and GU) is not enough to meet the dispatched power, the VPP will export as much power as resource availability allows.

    • The VPP can enter an idle state under three conditions: (ⅰ) the dispatched active power is equal to zero; (ⅱ) the available power within the VPP is zero, which would make the VPP unable to meet any dispatched power, and (ⅲ) SOUT (the apparent power at the AC-side) is smaller than the inverter minimum power. In idle state the VPP will only draw active power, determined by the user-defined No-load parameter and ZIP coefficients.

    • The last step in the procedure verifies that SOUT is not greater than the inverter rated power; if so, the model will adjust POUT, QOUT, PDG, and PES so the resulting SOUT is equal to the inverter rated power. This function is also used under the Follow operation mode.

    • For a time-driven simulation, the procedure must be repeated for every time-step.


    3. Simulation of the VPP Model Under Dispatch Operation Mode


    3.1. Test system configuration and characteristics

    This section is aimed at illustrating the behavior of the VPP model running under Dispatch operation mode. Figure 3 shows the test system: it is a three-phase 60-Hz overhead system serving two spot loads supplied from a HV/MV substation transformer. The VPP is connected to the MV network through a distribution MV/LV transformer.

    Figure 3. Case study 1: Test system configuration.

    The main characteristics of the test system are as follows:

    • HV rated voltage: 230 kV

    • MV rated voltage: 4.16 kV

    • LV rated voltage: 0.4 kV

    • Substation transformer rated power: 500 kVA

    • MV/LV transformer rated power: 210 kVA.

    System loads have been modeled as voltage-independent. Load profiles, solar irradiance, wind speed, and panel and ambient temperature curves were derived with the procedure presented in [14]. The actual voltage at the substation HV terminals is 1.05 pu.

    The main characteristics of the VPP are presented in Table 2; the study will consider both the Solar and Wind VPP. The efficiency curves used for the inverter and ES unit are shown in Figures 4 and 5, the power vs. Temperature curve of the PV generator is presented in Figure 6, and the power vs. wind speed curve is shown in Figure 7.

    Table 2. Case study 1: VPP characteristics.
    PV Generator Energy Storage Unit
    Rated power: 200 kW Rated power: 200 kW
    Minimum power: 5% Rated energy: 2000 kWh
    Power curtailment: 100% Initial stored energy: 1000 kWh
    Minimum stored energy: 400 kWh
    Minimum power: 5%
    Wind Generator Inverter
    Rated power: 200 kW Rated power: 200 kVA
    Hub height: 35 m Rated voltage: 0.4 kV
    Cut-in speed: 4 m/s Power factor: 0.9
    Cut-out speed: 20 m/s Minimum power: 5%
    Wind Reference height: 50 m No-load power: 1%
    Altitude over sea-level: 64 m Number of phases: 3
    Roughness length in wind direction: 0.8 m Connection: delta
    Minimum power: 2.4% ZIP coefficients: [0.9, 0, 0.1]
    Power curtailment: 100%
     | Show Table
    DownLoad: CSV
    Figure 4. Case study 1: Inverter efficiency curve.
    Figure 5. Case study 1: ES unit efficiency curve.
    Figure 6. Case study 1: Power vs. temperature curve.
    Figure 7. Case study 1: Power vs. wind speed curve.

    3.2. Simulation results

    The test system has been simulated for a period of 24 hours with a 1-hour time-step and system curves generated for a 1-year period. The selected 24 hours belong to a specific day of the year. Five different scenarios have been considered for simulation:

    (1) system without VPP (i.e. without DG and ES), also referred to as base case;

    (2) system with PV generation only;

    (3) system with Solar VPP;

    (4) system with only wind generation;

    (5) system with Wind VPP.

    The day under consideration presents distinctive variations between high and low load conditions (see Figure 8a); moreover, there is no significant correlation between system load and renewable resources, which will help to accentuate the advantages of a VPP with respect to DG without ES. Note that the load curves for the previous and following day present different values and characteristics.

    Figure 8. Case study 1: Base case.

    Figure 8b depicts the load voltages for the base case. For Scenarios 3 and 5, the VPPs operate under the Dispatch operation mode and follow the dispatch curve depicted in Figure 9. Take into account that the dispatch curve has not been optimized and no reactive power is dispatched in these scenarios. The operating conditions resulting from all scenarios are shown in Table 3.

    Figure 9. Case study 1: VPP dispatch curve.
    Table 3. Case study 1: Test system operating conditions.
    Base case Only PV Solar VPP Only Wind Wind VPP
    Active energy served (kWh) 1 4206.1 2851.7 2747.4 2852.9 2661.1
    Reverse energy flow (kWh) 1 - –104.4 0.0 –64.4 0.0
    Energy losses (kWh) 2 76.5 80.5 62.0 76.2 60.9
    Maximum active power (kW)1 238.1 226.9 176.7 213.9 145.8
    Minimum active power (kW)1 121.1 –43.8 75.9 –49.1 75.9
    Maximum voltage (pu) 1.0500 1.0556 1.0500 1.0593 1.0500
    Minimum voltage (pu) 0.9902 0.9971 1.0017 0.9953 1.0083
    DG/VPP energy (kWh)3 - 1462.9 1444.3 1417.2 1529.4
    1 Measured at the MV-side of the substation transformer.
    2 Measured losses do not include losses in the substation transformer.
    3 Measured at the LV-side of the interconnection transformer.
     | Show Table
    DownLoad: CSV

    Figure 10 depicts the active power measured at the MV-side of the substation transformer. As expected, there is a reduction in the active energy supplied by the substation transformer in Scenarios 2 through 5; moreover, all these scenarios also produced improvements in peak power and minimum voltage (see Figures 10 and 11). This is due to the fact that the peak load occurs in the morning when PV and wind generation are non-zero; however, there is a second peak at nighttime, on which DG without ES has little or no effect.

    Figure 10. Case study 1: Active power curves measured at MV-side substation terminals.
    Figure 11. Case study 1: Load node voltages-Phase A.

    Another important aspect is the effect of the intermittent nature of the solar and wind generation: in Scenarios 2 and 4 DG is allowed to produce as much energy as possible, but due to the low correlation between generation and load, too high system voltages are produced.

    The VPPs were capable of producing a greater impact on the operating conditions as both can adjust the power output according to the load behavior. The resulting active power and load voltage profiles present a much flatter behavior than the other scenarios (see Figures 10 and 11).

    The VPP dispatch curve was designed to mirror the behavior of the active power curve in the base case; consequently, the VPP power output is a reflection of the system's state and undesired conditions are avoided (e.g. too high system voltages).

    For the present case study, the Wind VPP was capable of following the dispatch curve without any deviations, whereas the Solar VPP failed to meet the dispatched power at the 9th and 10th hours, see Figure 12. The solar irradiance curve and the demanded power in the early hours of the day, forced the Solar VPP to depend entirely on the energy stored in the ES unit; by the 8th hour the stored energy had reached its minimum value and the power produced by the PV generator was not enough to meet the dispatch curve. As of the 11th hour, the power produced by the PV generator was sufficient to satisfy the dispatched power and charge the ES unit for the upcoming night hours; the operation curves of the VPP components are presented in Figure 13.

    Figure 12. Case study 1: Active power output.
    Figure 13. Case study 1: VPP components operation curves.

    On the other hand, the wind profile allowed the Wind VPP to charge the ES unit with enough energy to meet the prescheduled power output for the entire day. Note that the wind generation in the second half of the day is much lower when compared to that in the first half; however, the proper operation of the VPP permitted to store this excess of energy so it could be released when it was required later on.

    The VPP operates under the idle state in two different periods, as determined by the dispatch curve. Information related to these time periods is presented in Table 4; it can be observed that the power absorbed by the VPP is not constant due to the voltage dependency of the VPP ZIP load model. It is worth noting that when GU production is non-zero, the entirety of this power is used to charge the ES unit. Another important aspect is that a power curtailment is carried out in the Wind VPP at the 5th hour; due to the different corrections, the wind generated power is slightly greater than the GU rated power (201.3 kW). However, the ES unit can absorb a maximum power of 200 kW; as a consequence, the wind power was curtailed to meet the ES unit requirements.

    Table 4. Case study 1: VPP idle state.
    Solar VPP Wind VPP
    Hour Dispatched POUT Idle POUT PDG PES Idle POUT PDG PES
    4 0.0 2.0913 0.0 0.0 2.0913 103.4 –103.4
    5 0.0 2.0855 0.0 0.0 2.0853 200.0 –200.0
    16 0.0 2.0867 163.6 –163.6 2.0867 0.0 0.0
    17 0.0 2.0968 127.4 –127.4 2.0968 0.0 0.0
    18 0.0 2.0962 67.3 –67.3 2.0963 42.1 –42.1
    19 0.0 2.0933 21.5 –21.5 2.0933 31.9 –31.9
     | Show Table
    DownLoad: CSV

    4. Day-Ahead Dispatch Optimization


    4.1. Optimization procedure

    The dispatch curve used in the previous section was not derived by means of any optimization procedure; however, in order to maximize the potential benefits to the distribution system it is necessary for the VPP to follow a curve that has been optimized for the evaluation period under consideration. The VPP operation can be conducted from the perspective of a utility (i.e. aiming to improve network conditions) or operated by a private operator (i.e. with the objective of generating profits). The VPP dispatch optimization can be implemented following an approach similar to those used for determining the optimum operation of ES (e.g. [18,19,20,21,22,23,24,25,26,27]).

    Day-ahead dispatch is based on the system's forecast for the following day; the VPP operation must be optimized according to these forecasted load and renewable resource curves. This section presents a procedure for determining the optimal hourly day-ahead dispatch of multiple VPPs from the utility's perspective.

    In this section the VPP optimal dispatch is determined by means of a general parallel genetic algorithm (GPGA). The implemented GPGA is a global single-population GA [28], where only one population is created and the evaluation of each member of a given generation is independent from the rest of the members; therefore, the evaluation of a generation can be seen as an embarrassingly parallel problem [29]. This approach yields a reduction of the total execution times by distributing the evaluation of individuals in a generation among available cores in a multi-core installation.

    The GPGA will seek to determine the optimal values for every point of the VPP operation curve that minimize the following operation variables:

    • Distribution energy losses, without including losses in substation transformer (EL).

    • Energy supplied from the MV-side of the substation transformer (E).

    • Peak power supplied by the substation transformer from the MV-side terminals (PMAX).

    For a day-ahead dispatch with a 1-hour time step, the VPP operation curve will consist of 24 values, which can be either positive or negative (according to the power injection or absorption). Note that the procedure can cope with multiple dispatch curves; that is, different VPPs can follow different operation curves. The parameters used in the optimization procedure are shown in Table 5.

    Table 5. Case study 2: GPGA parameters.
    Maximum number of generations 200
    Number of traits 24 × Number of dispatch curves
    Population size 10 × Number of traits
    Epsilon (ε)-Tolerance 0.0001
    Delta (δ)-Number of generations for best fitness checking 20
    Crossover probability (pc) 0.80
    Mutation probability (pm) 0.02
    Elitism Yes
    Maximum curve value 1 1.00
    Minimum curve value1 –1.00
    1 It is a factor of the VPP inverter rated power.
     | Show Table
    DownLoad: CSV

    Two stopping criteria have been established for the optimization procedure:

    (1) The generation count (i.e. the number of executed generations) is equal to the established maximum number of generations.

    (2) The best fitness changes less than ε over δ generations.

    The GPGA has been implemented in MATLAB and the evaluation of each member in a population is conducted by simulating the test system during one day using a 1-hour time-step. Individual executions are distributed among 4 cores using the library developed by M. Buehren [30]. The procedure was tested on a laptop computer with an Intel Core i7-6700HQ processor (4 Cores, Clock frequency = 2.6–3.5 GHz), 12 GB RAM, and Windows 10 OS.

    The fitness function calculated for every evaluation will be equal to the square root of the sum obtained from the operation variables:

    Fitness=(EVPPEBASE)2+(ELVPPELBASE)2+(PMAX_VPPPMAX_BASE)2+(ΔRPF)2+NVPPi=1ΔE2i+NVPPi=1CD2i (8)
    ΔRPF={0,PMIN_VPP>0PMIN_VPPPMAX_BASE,PMIN_VPP0 (9)
    ΔE={0,EES_LIMEESEES_LIM+EES_OBJEESEES_OBJ,else (10)

    where NVPP is the number of VPPs in the system, PMIN is the minimum power served by the substation transformer from the MV-side; EES_OBJ is the expected stored energy in the ES unit at the end of the evaluation period; EES_LIM and EES_LIM+ are the lower and upper limits of a bandwidth around EES_OBJ; CD is the Euclidean distance between the proposed dispatch curve and the actual VPP output curve. The sub index BASE makes reference to the base case (i.e. system without DG or VPP), while VPP is used for those quantities obtained when the VPPs are present in the system. ∆RPF has been introduced in the fitness function as penalty for those cases in which reverse power flow occurs, while ∆E is a penalty for those cases when the ES unit stored energy falls outside of the bandwidth defined by EES_LIM and EES_LIM+. Forcing the ES unit to maintain a minimum value of stored energy at the end of the evaluation period will ensure that the VPP will have enough energy for meeting the following day's dispatch (i.e. the day after the day-ahead).

    The procedure seeks to minimize the value of the fitness function presented in Equation (8). For every evaluation, the procedure checks that all voltages are within accepted limits (between 0.95 and 1.05 pu), no network element experiences overload, and the iterative algorithm converges for all solution steps. If one or more of these conditions are not fulfilled, the procedure will assign a large value to the fitness function, which will make difficult its participation in the subsequent generations of the GPGA.


    4.2. Test system configuration

    The new test system prepared for this case study is a three-phase 60-Hz overhead system with a simplified representation of the HV transmission system (see Figure 14). VPPs use the same configuration as in the previous test system, see Figures 3 through 7. By default, all distribution loads are connected to the system through MV/LV transformers.

    Figure 14. Case study 2: Test system configuration.

    Some important information related to the new test system follows:

    • HV rated voltage: 230 kV.

    • HV actual voltage: 1.03 pu.

    • MV rated voltage: 4.16 kV.

    • LV rated voltage: 0.4 kV.

    • Substation transformer rated power: 2500 kVA.

    • Total rated active load: 1700 kW.

    As in Test System 1, loads have been modelled as voltage-independent and curve shapes have been derived with the procedure presented in [14].

    Solar VPP 1 and Wind VPP 1 use the same parameters as those in case study 1; furthermore, the expected stored energy, and lower and upper bandwidth limits have been set to 1500,800, and 1200 kWh, respectively. The characteristics of Solar VPP 2 and Wind VPP 2 are presented in Table 6 (all VPPs operate under the Dispatch operation mode). Note that the rated values and locations of all VPPs have been arbitrarily chosen. It is assumed that the solar irradiance and reference wind speed are the same for the entire system; due to height corrections the actual wind speed curve will be different for both Wind VPPs.

    Table 6. Case study 2: VPP characteristics.
    Solar VPP 2
    PV Generator Energy Storage Unit Inverter
    Rated power: 300 kW Rated power: 300 kW Rated power: 300 kVA
    Minimum power: 5% Rated energy: 3000 kWh Rated voltage: 0.4 kV
    Power curtailment: 100% Initial stored energy: 1500 kWh Power factor: 0.9
    Minimum stored energy: 600 kWh Minimum power: 5%
    Minimum power: 5% No-load power: 1%
    Expected stored energy: 1500 kWh Number of phases: 3
    Lower limit: 1200 kWh Connection: delta
    Upper limit: 1800 kWh ZIP coefficients: [0.9, 0, 0.1]
    Wind VPP 2
    Wind Generator Energy Storage Unit Inverter
    Rated power: 300 kW Rated power: 300 kW Rated power: 300 kVA
    Hub height: 45 m Rated energy: 3000 kWh Rated voltage: 0.4 kV
    Cut-in speed: 4 m/s Initial stored energy: 1500 kWh Power factor: 0.9
    Cut-out speed: 20 m/s Minimum stored energy: 600 kWh Minimum power: 5%
    Wind Reference height: 50 m Minimum power: 5% No-load power: 1%
    Altitude over sea-level: 64 m Expected stored energy: 1500 kWh Number of phases: 3
    Roughness length: 0.8 m Lower limit: 1200 kWh Connection: delta
    Minimum power: 2.4% Upper limit: 1800 kWh ZIP coefficients: [0.9, 0, 0.1]
    Power curtailment: 100%
     | Show Table
    DownLoad: CSV

    4.3. Simulation results

    The GPGA was executed to obtain the optimal day-ahead dispatch of active power for the VPPs in the system. Due to the different nature of wind and solar resources, the procedure will derive two dispatch curves: the first one will be assigned to the Solar VPPs, while the Wind VPPs will operate according to the second dispatch curve. As in the previous study, the selected 24 hours belong to a specific day of the year, and the load curves for the previous and following day present different values and characteristics.

    A GA is a stochastic optimization procedure; thus, it is possible for it to converge before reaching the optimum solution; therefore, two different scenarios have been considered for the GPGA execution:

    (1) the GPGA will consider both stopping criteria mentioned above;

    (2) the procedure is executed without the δ/ε criterion; that is, the procedure will stop when it has reached the maximum number of generations.

    The results for both scenarios are presented in Table 7. The table also includes the operation variables for the base case. These results show a clear improvement on the energy and peak power served by the substation transformer; both operation variables (as well as the best fitness) present consistent values for both scenarios. On the other hand, it was not possible to produce a reduction on the energy losses; however, the energy losses for the base case do not include the losses of the VPP interconnection transformers, which represent a new and significant source of losses; energy losses including only the interconnection transformers amount to 618.4 kWh.

    Table 7. Case study 2: Simulation results.
    Base case Scenario 1 Scenario 2
    Best fitness - 1.3827 1.3765
    Generation count - 57 200
    Simulation time (sec) - 501.7 1625.0
    Energy (kWh) 22021.2 12722.6 12753.5
    Energy losses (kWh) 566.6 583.4 589.4
    Maximum power (kW) 1178.9 792.4 752.7
    Stored energy Solar VPP 1 (kWh)1 - 828.4 808.4
    Stored energy Solar VPP 2 (kWh)1 - 1242.6 1212.6
    Stored energy Wind VPP 1 (kWh)1 - 744.7 764.7
    Stored energy Wind VPP 2 (kWh)1 - 1763.9 1779.9
    1 Measured at the end of the evaluation period.
     | Show Table
    DownLoad: CSV

    Figure 15 depicts the active power profile for the base case and the resulting profile from Scenario 1; additionally, it presents the obtained dispatch curves and the actual VPP operation curve. Figure 15a shows that the VPP operation produced a sustained reduction on the active power served from the substation transformer, while Figures 15b and 15c demonstrate that the VPPs are capable of accurately following their dispatch curves.

    Figure 15. Case study 2: Simulation results.

    The only deviation present in the VPP dispatch occurs for the Solar VPP; the dispatch curve set a charging power of –0.02 pu at the 22nd hour; however, this value is smaller than the inverter minimum power, which forces the VPP to enter the idle state. It is important to observe that both dispatch curves present different behaviors throughout the evaluation period, which is likely to be caused by the different nature of both resources.

    Due to difference in hub height, the Wind VPP 2 will produce more power and energy (in pu) than the Wind VPP 1; as a result, the stored energy in the ES unit at the end of the evaluation period (with respect to the rated energy) will also be greater, see Table 7. In general, the stored energy in the Wind VPP 2 at the end of the evaluation period is always close to the bandwidth's upper limit, which means that the VPP could inject more energy into the system and still remain within the desired bandwidth; however, the Wind VPP 2 cannot use this extra energy if it is forced to follow the same dispatch curve as the Wind VPP 1. In order to explore the possibility of utilizing the extra energy stored in the Wind VPP 2, the GPGA has been executed considering one dispatch curve for each Wind VPP. Additionally, the GPGA was also executed considering one dispatch curve per VPP (i.e. 4 dispatch curves). Simulation results for these new conditions are presented in Table 8.

    Table 8. Case study 2: Simulation results-Multiple dispatch curves.
    3 Dispatch curves1 4 Dispatch curves2
    Scenario 1 Scenario 2 Scenario 1 Scenario 2
    Best fitness 1.3646 1.3233 1.3936 1.3999
    Generation count 52 200 23 200
    Simulation time (sec) 680.8 2547.9 364.9 3382.1
    Energy (kWh) 12696.1 12461.8 12941.9 12800.4
    Energy losses (kWh) 596.7 584.5 596.2 622.4
    Maximum power (kW) 758.9 711.3 803.8 757.1
    Stored energy Solar VPP 1 (kWh)1 860.4 827.8 802.3 840.7
    Stored energy Solar VPP 2 (kWh)1 1290.5 1241.7 1762.4 1211.5
    Stored energy Wind VPP 1 (kWh)1 832.3 890.2 964.2 860.9
    Stored energy Wind VPP 2 (kWh)1 1425.8 1288.7 1215.3 1301.9
    1 One dispatch curve per Wind VPP.
    2 One dispatch curve per VPP (Solar and Wind).
     | Show Table
    DownLoad: CSV

    Results in Table 8 show that by defining a dispatch curve for each Wind VPP (i.e. 3 dispatch curves) a slight improvement in the energy and maximum power served by the substation transformer has been produced; moreover, it allowed the Wind VPP 2 to inject as much energy as possible into the system (while respecting the stored energy bandwidth). On the other hand, using one dispatch curve per VPP (i.e. 4 dispatch curves) produced no improvement with respect to the previous executions of the optimization procedure.

    Tables 7 and 8 show that Scenario 2 produced better results than Scenario 1 (in terms of best fitness values), although the operation variables present low variations; moreover, using 3 dispatch curves produced the best results among all executions. This behavior indicates that VPPs with different characteristics and GU output (in pu) can be operated using individual dispatch curves.

    Simulation times have been affordable for all GPGA executions, ranging from 5 minutes to 1 hour; these times vary according to the number of executed generations and population size, which is a function of the number of dispatch curves to be optimized (see Table 5).


    5. Conclusions

    This paper has presented a virtual power plant model for power flow calculations. The model has been compiled as a custom-made capability within the OpenDSS COM DLL and can be used in time-driven simulations controlled from other software platforms, e.g. MATLAB. The model has the same capabilities as other OpenDSS objects and it can be defined, modified, and accessed using the same methods.

    Two different operation modes have been implemented inside the model: Follow and Dispatch. Under the Follow mode the model behaves according to the curves assigned to the GU and ES unit, whereas in Dispatch mode the VPP acts as a dispatchable unit.

    The list of features implemented in the model includes: (ⅰ) active and reactive power dispatch, (ⅱ) power curtailment, (ⅲ) inverter idle power modeled as a ZIP load, (ⅳ) minimum power limit for GU and ES unit, (ⅴ) automatic calculation of ES unit maximum charge/discharge power, and (ⅵ) inverter power limits.

    The first case study has shown how the VPP can make better use of the renewable resources, since it is capable of storing the generation surplus and injecting it into the network during low generation periods; the presence of ES is essential for meeting the power dispatch curve.

    Contrary to DG based only on intermittent resources, the VPP output can be adjusted to meet the system requirements and thus avoid undesired operating conditions, such as system voltages above accepted limits. This is a vital feature for achieving larger penetration factors of renewable generation without compromising normal system operation.

    The VPP can help to further improve system operating conditions by injecting or absorbing reactive power; the absorption of reactive power is expected to play a key role if the objective is to export energy to the transmission system [31].

    The paper has also presented a procedure for determining the optimal day-ahead dispatch for the VPP based on a GPGA. The procedure seeks to explore the potential improvement on the overall operating conditions when the VPP is operated from the perspective of a utility. The presented methodology relies on a power flow simulator (i.e. OpenDSS) for calculating the system operating conditions; namely, no simplifications are introduced into the solution of the electrical system. Furthermore, the evaluated system can be of any size and include detailed models and control algorithms for system components (e.g. the presented VPP model). An accurate calculation of system conditions is essential for the optimization methodology since a poor estimation would cause the procedure to produce erroneous results.

    A genetic algorithm requires a large number of executions in order to produce accurate results; although each individual execution does not require large simulation times, the repeated evaluation of the system could result in prohibitive executions times if they were to be performed in a sequential manner (i.e. using single core computing). This work makes use of parallel computing in combination with the GA (i.e. a GPGA) in order to distribute the individual executions among the available cores to produce a reduction in the total execution times.

    The results presented in the second case study show that the proposed optimization procedure is capable of producing significant improvements on the system conditions and that the parallel approach used in this work produces affordable simulation times. These simulation times could be further reduced through the use of a larger multi-core installation.

    Similar results to those presented in this paper could be achieved by defining appropriate operation curves of an OpenDSS Storage object that works in conjunction with the PV System or Generator objects; however, the coordination between both objects should be handled by an external tool (e.g. MATLAB), which would represent a significant loss in performance. The main advantage of the present model is that all operations and control decisions are compiled within the COM DLL, which makes them seamless to the user and helps to maintain the high computing performance of OpenDSS.

    In its present form the developed VPP model is adequate for exploring the effect that the aggregation of renewable generation and energy storage can have on the system. Future work could be aimed at expanding the VPP model to incorporate losses in the DC network and implement an algorithm to automatically control the inverter's reactive power based on the system's conditions. Furthermore, an optimization procedure must be implemented in order to determine the optimal VPP operation curve for an extended evaluation period (i.e. one year or more).


    Conflict of Interest

    The author declares no conflicts of interest in this paper.




    [1] S. Datta and D. Nettleton, Statistical Analysis of Next Generation Sequencing Data, Springer, 2014. Available from: http://link.springer.com/content/pdf/10.1007/978-3-319-07212-8.pdf.
    [2] W. Sun and Y. Hu, Mapping of expression quantitative trait loci using RNA-seq data, in Statistical Analysis of Next Generation Sequencing Data (eds. D. Nettleton and S. Datta), 2014, 25–50.
    [3] P. S. Schnable and N. M. Springer, Progress toward understanding heterosis in crop plants, Annu. Rev. Plant Biol., 64 (2013), 71–88.
    [4] A. Paschold, Y. Jia, C. Marcon, et al., Complementation contributes to transcriptome complexity in maize (Zea mays L.) hybrids relative to their inbred parents., Genome Res., 22 (2012), 2445–2454.
    [5] G. D. M. Bell, N. C. Kane, L. H. Rieseberg, et al., RNA-Seq analysis of allele-specific expression, hybrid effects, and regulatory divergence in hybrids compared with their parents from natural populations, Genome Biol. Evol., 5 (2013), 1309–1323.
    [6] J. K. Pickrell, J. C. Marioni, A. A. Pai, et al., Understanding mechanisms underlying human gene expression variation with rna sequencing, Nature, 464 (2010), 768–772.
    [7] W. Sun and Y. Hu, eQTL Mapping Using RNA-seq Data, Stat. Biosci., 5 (2013), 198–219.
    [8] C. T. Harvey, G. A. Moyerbrailean, G. O. Davis, et al., Quasar: quantitative allele-specific analysis of reads, Bioinformatics, 31 (2014), 1235–1242.
    [9] N. Raghupathy, K. Choi, M. J. Vincent, et al., Hierarchical analysis of RNA-seq reads improves the accuracy of allele-specific expression, Bioinformatics, 34 (2018), 2177–2184.
    [10] S. Srivastava and L. Chen, A two-parameter generalized Poisson model to improve the analysis of RNA-seq data., Nucleic Acids Res., 38 (2010), e170.
    [11] X. Wei and X. Wang, A computational workflow to identify allele-specific expression and epigenetic modification in maize., Genom. Proteom. Bioinf., 11 (2013), 247–252.
    [12] M. D. Robinson, D. J. McCarthy and G. K. Smyth, edgeR: a Bioconductor package for differential expression analysis of digital gene expression data., Bioinformatics (Oxford, England), 26 (2010), 139–140.
    [13] D. J. Lorenz, R. S. Gill, R. Mitra, et al., Using RNA-seq Data to Detect Differentially Expressed Genes, in Statistical Analysis of Next Generation Sequencing Data (eds. S. Datta and D. Nettleton), 2014, chapter 2, 25–49.
    [14] Y.-J. Hu, W. Sun, J.-Y. Tzeng, et al., Proper use of allele-specific expression improves statistical power for cis -eQTL mapping with RNA-seq data, J. Am. Stat. Assoc., 110 (2015), 962–974.
    [15] W. Landau, J. Niemi and D. Nettleton, Fully bayesian analysis of rna-seq counts for the detection of gene expression heterosis, J. Am. Stat. Assoc., 114 (2019), 601–612.
    [16] N. I. Panousis, M. Gutierrez-Arcelus, E. T. Dermitzakis, et al., Allelic mapping bias in RNA-sequencing is not a major confounder in eQTL studies, Genome. Biol., 15 (2014), 467.
    [17] J. F. Degner, J. C. Marioni, A. A. Pai, et al., Effect of read-mapping biases on detecting allele-specific expression from RNA-sequencing data, Bioinformatics, 25 (2009), 3207–3212.
    [18] R. Vijaya Satya, N. Zavaljevski and J. Reifman, A new strategy to reduce allelic bias in RNA-Seq readmapping, Nucleic Acids Res., 40 (2012), 1–9.
    [19] K. R. Stevenson, J. D. Coolon and P. J. Wittkopp, Sources of bias in measures of allele-specific expression derived from RNA-sequence data aligned to a single reference genome., BMC Genom., 14 (2013), 536.
    [20] Y. Chen, A. T. L. Lun and G. K. Smyth, Differential expression analysis of complex RNA-seq experiments using edgeR, in Statistical Analysis of Next Generation Sequencing Data, Springer, Cham, 2014, 51–74.
    [21] T. Park and G. Casella, The Bayesian lasso, J. Am. Stat. Assoc., 103 (2008), 681–686.
    [22] C. M. Carvalho, N. G. Polson and J. G. Scott, Handling Sparsity via the Horseshoe, J. Mach. Learn. Res., 5 (2009), 73–80.
    [23] A. Gelman, Prior distributions for variance parameters in hierarchical models, Bayesian Analysis, 1 (2006), 515–533.
    [24] A. Gelman, J. B. Carlin, H. S. Stern,et al., Bayesian Data Analysis, CRC press, 2013.
    [25] J. K. Ghosh, M. Delampady and T. Samanta, An Introduction to Bayesian Analysis, Springer, 2006. Available from: http://onlinelibrary.wiley.com/doi/10.1002/9781118684818.ch16/summary.
    [26] L. G. León-Novelo, L. M. McIntyre, J. M. Fear, et al., A flexible Bayesian method for detecting allelic imbalance in RNA-seq data, BMC Genom., 15 (2014), 920.
    [27] J. Niemi, E. Mittman, W. Landau, et al., Empirical Bayes analysis of RNA-seq data for detection of gene expression heterosis, J. Agr. Biol. Envir. St., 20 (2015), 614–628.
    [28] M. A. Van De Wiel, G. G. R. Leday, L. Pardo, et al., Bayesian analysis of RNA sequencing data by estimating multiple shrinkage priors, Biostatistics, 14 (2013), 113–128.
    [29] W. Landau and J. Niemi, A fully Bayesian strategy for high-dimensional hierarchical modeling using massively parallel computing, 2016. Available from: http://arxiv.org/abs/1606.06659.
    [30] A. Lithio and D. Nettleton, Hierarchical modeling and differential expression analysis for RNA-seq experiments with inbred and hybrid genotypes, J. Agr. Biol. Envir. St., 20 (2015), 598–613.
    [31] M. Ventrucci, E. M. Scott and D. Cocchi, Multiple testing on standardized mortality ratios: a Bayesian hierarchical model for FDR estimation, Biostatistics, 12 (2011), 51–67.
    [32] P. Muller, G. Parmigiani and K. Rice, FDR and Bayesian multiple comparisons ules, 2006. Available from: http://biostats.bepress.com/jhubiostat/paper115.
    [33] H. Y. Bar, J. G. Booth and M. T. Wells, A bivariate model for simultaneous testing in bioinformatics data, J. Am. Stat. Assoc., 109 (2014), 537–547.
    [34] P. Müller, G. Parmigiani, C. Robert, et al., Optimal sample size for multiple testing: the case of gene expression microarrays, J. Am. Stat. Assoc., 99 (2004), 990–1001.
    [35] S. Anders and W. Huber, Differential expression analysis for sequence count data, Genome Biol., 11 (2010), R106.
    [36] P. R. Hahn and J. He, Elliptical slice sampling for Bayesian shrinkage regression with applications to causal inference, 2016. Available from: http://faculty.chicagobooth.edu/richard. hahn/JCGS_submit.pdf.
    [37] M. C. Sachs, plotROC: A tool for plotting roc curves, J. Stat. Software, 79 (2017), 1–19.
    [38] Y. Benjamini and Y. Hochberg, Controlling the false discovery rate: A practical and powerful approach to multiple testing, J. R. Stat. Soc-B, 57 (1995), 289–300.
    [39] P. S. Schnable, D. Ware, R. S. Fulton, et al., The B73 maize genome: complexity, diversity, and dynamics, Science, 326 (2009), 1112–1115.
  • This article has been cited by:

    1. Talal Alzahrani, Raluca Eftimie, Dumitru Trucu, Multiscale moving boundary modelling of cancer interactions with a fusogenic oncolytic virus: The impact of syncytia dynamics, 2020, 323, 00255564, 108296, 10.1016/j.mbs.2019.108296
    2. MONIKA JOANNA PIOTROWSKA, MAREK BODNAR, URSZULA FORYŚ, Tractable Model of Malignant Gliomas Immunotherapy with Discrete Time Delays, 2014, 21, 0889-8480, 127, 10.1080/08898480.2013.804690
    3. Alberto Gandolfi, Andrea Pugliese, Carmela Sinisgalli, Epidemic dynamics and host immune response: a nested approach, 2015, 70, 0303-6812, 399, 10.1007/s00285-014-0769-8
    4. Monika Joanna Piotrowska, Marek Bodnar, Jan Poleszczuk, Urszula Foryś, Mathematical modelling of immune reaction against gliomas: Sensitivity analysis and influence of delays, 2013, 14, 14681218, 1601, 10.1016/j.nonrwa.2012.10.020
    5. Marcello Delitala, Umberto Dianzani, Tommaso Lorenzi, Matteo Melensi, A mathematical model for immune and autoimmune response mediated byT-cells, 2013, 66, 08981221, 1010, 10.1016/j.camwa.2013.06.026
    6. Todd R. Young, Richard Buckalew, Addison K. May, Erik M. Boczko, A low dimensional dynamical model of the initial pulmonary innate response to infection, 2012, 235, 00255564, 189, 10.1016/j.mbs.2011.12.004
    7. Todd R Young, Erik M Boczko, Early treatment gains for antibiotic administration and within human host time series data, 2018, 35, 1477-8599, 203, 10.1093/imammb/dqw025
    8. Monika J. Piotrowska, Urszula Foryś, The nature of Hopf bifurcation for the Gompertz model with delays, 2011, 54, 08957177, 2183, 10.1016/j.mcm.2011.05.027
    9. Shujing Shi, Jicai Huang, Jing Wen, Shigui Ruan, Bifurcation Analysis of a Dynamical Model for the Innate Immune Response to Initial Pulmonary Infections, 2020, 30, 0218-1274, 2050252, 10.1142/S0218127420502521
    10. Paul A. Valle, Luis N. Coria, Diana Gamboa, Corina Plata, Bounding the Dynamics of a Chaotic-Cancer Mathematical Model, 2018, 2018, 1024-123X, 1, 10.1155/2018/9787015
    11. Florent Feudjio Kemwoue, Vandi Deli, Joseph Marie Mendimi, Carlos Lawrence Gninzanlong, Jules Fossi Tagne, Jacques Atangana, Dynamics of cancerous tumors under the effect of delayed information: mathematical and electronic study, 2022, 2195-268X, 10.1007/s40435-022-01031-2
    12. Florent Feudjio Kemwoue, Vandi Deli, Hélène Carole Edima, Joseph Marie Mendimi, Carlos Lawrence Gninzanlong, Mireille Mbou Dedzo, Jules Fossi Tagne, Jacques Atangana, Effects of delay in a biological environment subject to tumor dynamics, 2022, 158, 09600779, 112022, 10.1016/j.chaos.2022.112022
  • Reader Comments
  • © 2019 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5901) PDF downloads(465) Cited by(1)

Article outline

Figures and Tables

Figures(10)  /  Tables(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog