# Microtubule Motility Analysis based on Time-Lapse Fluorescence Microscopy

S. Masoudia , C. H.G. Wrighta , J. C. Gatlinb and J. S. Oakeyc
a Electrical and Computer Engineering Department
University of Wyoming, Laramie, WY 82071
b Molecular Biology Department
University of Wyoming, Laramie, WY 82071
c Chemical Engineering Department
University of Wyoming, Laramie, WY 82071

###### Abstract

This paper describes an investigation into part of the mechanical mechanisms
underlying the formation of mitotic spindle, the cellular machinery
responsible for chromosomal separation during cell division. In normal
eukaryotic cells, spindles are composed of microtubule filaments that radiate
outward from two centrosomes. In many transformed cells, however, centrosome
number is misregulated resulting in cells with more than two centrosomes.
Addressing the question of how these cells accommodate these additional
structures by coalescing supernumerary centrosomes to form normal spindles
will provide a powerful insight toward understanding the proliferation of
cancer cells and developing new therapeutics. The process of centrosome
coalescence is thought to involve motor proteins that function to slide
microtubules relative to one another. Here we use in vitro motility assays
combined with fluorescence microscopy to visualize, characterize and quantify
microtubule-microtubule interactions. After segmenting the microtubules, their
speed and direction of movement are the extracted features to cluster their
interaction type. In order to evaluate the potential of our processing
algorithm, we created a simulated dataset similar to the time-lapse series.
Once our procedure has been optimized using the simulated data, we will apply
it to the real data. Results of our analyses will provide a quantitative
description of interaction among microtubules. This is a potentially important
step toward more thorough understanding of cancer.

Keywords: Mitotic spindle, microtubule, centrosome, chromosome, fluorescence
microscopy, interaction, segmentation, tracking.

##  1 Introduction

Cancer is a somewhat generic word that refers to a category of all diseases
having out of control cell growth as their main characteristic. This non-stop
cell division leads to the formation of cell masses called tumors. Tumors can
interfere with the regular mechanism of digestive, nervous, and circulatory
systems or release hormones that change the body function [1, 2]. According to
the American Cancer Society, cancer is considered to be one of the leading
causes of morbidity and mortality worldwide [1]. It is the second most common
cause of death in the US, which accounts for nearly $1$ out of every $4$
deaths [1].

Studies over a century ago proposed aneuploidy (the presence of an abnormal
number of chromosomes in a cell) as the original source of cancer [3].
Aneuploidy is the most common genomic instability in human cancer cells [3].
The prevalence of centrosome amplification (i.e., in terms of number) in a
variety of tumor types and its direct correlation with aneuploidy, have
provided a strong foundation for at least one possible pathway to aneuploidy
in cancer cells [3].

To better understand the circumstances under which cancer cells might form,
one should first know about the mitotic spindle structure, components and
mechanism. Mitotic spindle is considered as a machine that segregates
chromosomes into two daughter cells during a certain type of cell division
called mitosis[3]. Main elements of the spindle are MT polymers, whose
intrinsic polarity and dynamic properties are crucial, for its proper function
[3]. In most cell types, MT nucleation occurs primarily at two spindle poles
called centrosomes [3, 4]. Normally, the opposing nature of dual centrosomes
facilitates the assembly of a bipolar spindle and ensures equal distribution
of the replicated genome to each daughter cell [3]. Under the abnormal
cancerous conditions, the overproduction of centrosomes leads to the formation
of multi-polar spindles, lagging chromosome and chromosomal instability (CIN)
[3, 5]. Many other factors, such as molecular motors and regulators of MT
dynamics, are associated with organizing the spindle [3].

The past decade has provided a wealth of information on the molecular organs
which are responsible for spindle assembly. Among this information, high-
resolution views of the detailed movements and dynamics of the spindle MTs and
chromosomes are of great interest [4]. A closer look into the mitotic spindle
revealed that MTs are its essential puzzle components. These polymers are
arranged head to tail to form a well-organized network of hollow tubes [6].
These tubes have two ends with different dynamic properties [6]. In the
cytoplasm environment, the minus ends are found throughout the spindle and
some are anchored to the centrosome while the more dynamic free plus ends
probe to reach specific targets during a search and capture process [6]. This
dynamic behavior of MTs is tightly regulated both spatially and temporally
[7]. The behavior itself, and its regulation are of great importance for
certain cellular functions such as cell division, signaling, polarity, shape
maintenance and so on [7]. Thus, modeling MTs’ dynamics and regulations is a
topic of major interest for researchers. Results of studying their dynamics,
movements, assembly, and how they interact with each other to form the
spindle, will potentially yield important insights into cancer cell biology
[8].

While significant investigations have been dedicated to this area, further
study on how exactly MTs participate in this process could shed more light on
new aspects of spindle assembly and, at a very basic level, the mechanistic
underpinnings of cancer. Former studies have attempted to track individual MTs
through in-vitro motility assays, disregarding the manner in which MT-MT
interaction influences the MT movements/dynamics. MTs must interact with each
other to fulfill several cellular functions, the most critical of which is
chromosome segregation during mitosis. A subset of MAPs have inherent ATPase
activity and function as motors which can walk along MTs. During mitotic
spindle assembly, two motors are particularly important: the minus-end
directed motor cytoplasmic dynein (dynein) and the plus-end directed kinesin
motor known as Eg-5. These motors bind to, move, and spatially organize MTs
relative to one another [9, 10, 11, 12]. Image analysis used in this study
should therefore inform us about MT velocities at each frame, elucidating the
function of MAPs and their contribution to MTs sliding, interaction, and
assembly [10, 11, 12].

For this purpose, $2$-D microscopy time-lapse images of MTs’ motility and
interaction were acquired. These images are obtained using time-lapse total
internal reflection fluorescence (TIRF) microscopy. The MTs in these images
are assembled in from Xenopus egg extracts spiked with fluorescently labeled
tubulin heterodimers [7] Currently, manual inspection of these time-lapse
images is the only available method, which is laborious, time consuming,
heavily biased and error-prone. Accordingly, there is a need for an automated
method which works much faster and provides better accuracy based on an
efficient low-cost algorithm. Image processing tools will greatly automate the
collection of a wide variety of metrics related to the MT-MT interactions.
Conventional image processing approaches in molecular cell biology, typically
consist of two subsequent steps of segmentation and object tracking.
Segmentation has been usually done based on thresholding, multiscale analysis
using the wavelet transform, or model fitting while nearest-neighbor or
smooth-motion criteria are basic methods for the second step of data
association [13]. Spatiotemporal-based methods and particle filtering have
also been proposed respectively for MT growth analysis which demonstrate to
work well [13]. Metrics of interest in this project include whether or not MTs
interact with each other while passing nearby. To accomplish this aim, two
main contributions are provided in this paper. First a new model for
simulating the MTs’ motility has been proposed that closely represents the
actual time-lapse images. Next, some standard deterministic methods have been
applied to this simulated data for performance evaluation of these algorithms
in a controlled way.

To substitute the manual inspection which is the currently used method for
detecting the potential MT-MT interaction in time-lapse images of MTs, an
automatic image processing algorithm is required. Image processing tools will
greatly automate the collection of a wide variety of metrics related to the
MT-MT interactions. Results of studying their dynamics, movements, assembly,
and how they interact with each other to form the spindle, will potentially
yield important insights into cancer cell biology.

The remainder of this paper is organized as follows. Section $2$ describes the
simulated data and introduces the processing algorithm. Section $3$ provides
the results generated by our processing algorithm on simulated data and
finally sections $4$ and $5$ present our conclusions.

##  2 Methods

###  2.1 Simulated Data

As was previously mentioned, the first step in this research project is to
create a simulated dataset which can approximately imitate the appearance and
dynamic behavior of MTs that are observed in the actual microscopy time-lapse
images. Similar to other cases where AI is applied to solve the biomedical
problems, main reason for generating the simulated data is to compensate for
the limited annotated data available for training. Besides, the simulated data
is used as input to the proposed image processing approaches which identify
individual MTs and track them in sequential time-lapse image frames. The
outcomes of this process will provide the basis with which to examine,
evaluate and make improvements to the processing algorithms.

Similar to the provided time-lapse series, the simulated dataset is a $2$-D+t
($x$ and $y$ spatial dimensions plus time) signal constructed with a sequence
of artificial images, each of which is assumed to be a sample at a certain
point in time. Considering the camera’s field of view to better resemble the
entrance and egress of the MTs to a time-lapse image frame, the central part
of a large predefined area is considered to be the main recorded frame. To
reproduce an environment with different conditions existing in the real time-
lapse images, the simulated data generator is coded in a flexible way so that
the number and size of the time-lapse image frames, time resolution, number of
participating MTs and the chance of MTs to interact are adjustable. To
assemble the final time-lapse image a specific number of produced frames are
attached sequentially.

Approximate simulation of the stochastic behavior of the MTs has been
fulfilled through several assumptions about their geometric specifications and
movement parameters. The geometric considerations mainly include the designed
structure of MTs. Each MT can be defined as a wagon train containing certain
number of wagons or blocks, where each block or wagon is formed by ${w\times
w}$ pixels. The first block is the head wagon while the last is its tail.
These wagons overlap to a certain amount that can be controlled by the code
generator. Other details are the length and width of an artificial MT. Overall
length is determined primarily by the number of wagons. For the length to be
realistically random, it is determined via a pseudo-random number generator
(PRNG) having a uniform probability distribution over an interval of
$(\frac{L}{2},L)$. $L$ limits the maximum possible length of the MTs.
Fortunately, all real-world MTs have basically the same width of $w$,
specified by the number of pixels along $x$ and $y$ axes in a block. Also
incorporated into the program is the real-world characteristic that length
might vary over time, which is achieved through defining $dL$ (the ratio of
maximum length change) and the uniform probability distribution over the
interval of $(-{dL},{dL})$.

Each MT has its own random initial conditions, including the coordinates of
the head block, initial direction and body shape. The initial location of the
head block follows a uniform probability distribution all around the main
frame. Direction only has four initial possible values:
${\\{{\frac{\pi}{4}},-{\frac{\pi}{4}},{\frac{3\pi}{4}},-{\frac{3\pi}{4}}\\}}$
radians.

Figure 1: A $288\times 512$ frame of (left) real frame obtained from sliding
MTs which belong to Xenopus egg extracts, (right) simulated data.

For reference, $0$ and $\pi/2$ radians would be seen as the left to right
horizontal and downward vertical directions, respectively. To obtain a random
initial direction, basic values are used to create the first frame. Then,
during the transition into the next few frames, stochastic moving parameters
help to turn MTs’ direction from a deterministic to a random nature.
Disregarding these first few frames for the final time-lapse image
compilation, results in random initial direction in the time-lapse image. This
would save us from adding to the number of possible initial directions and
employing a random initial body shape method at the same time. Using these
initial conditions, coordinates of each of the MT pixels are exactly
determined at the first frame.

Other aspects to the program are motility parameters that consist of its
direction and speed. Because real images of the MTs do not show them to have
highly “bent” shapes, available choices of moving directions should not be
more than $\pm{\frac{\pi}{4}}$. We implemented three allowed possibilities:
{straight, slightly right and slightly left}. Using a Gaussian probability
distribution with its mean value mapped to the straight direction, a random
direction of movement is obtained. Change of direction is allowed in every
five frames, again to improve the realism of the simulation. Speed also
follows a uniform probability distribution that is limited between user-
specified minimum and maximum values in units of pixels.

Finally, it should be noted that each pixel in a block that makes up that part
of the MT adds 32 levels of grayscale to the zero background (the full
grayscale range of the pixels is 256). Adding grayscales makes the inner parts
of the MT image lighter than its boundaries and also leads to much lighter
image regions at the areas over which MT-MT interactions occur. This contrast
provides greater fidelity to the real-world microscopy time-lapse images.

In our study, we generated 20 simulated sequences, each containing 379 frames
of size $256\times 256$ pixels. We used 4 randomly chosen sequences for the
test data and the rest for the training data.

###  2.2 Processing: segmentation

Image processing algorithms that are being used for this project can generally
be grouped into two categories: $1)$ those that process individual frames and
sequence them along the time axis later, or $2)$ those that process multiple
frames at the same time to track an object (a particular MT) along sequential
frames. The first method ignores any knowledge about possible correlation
among the following frames, and therefore may lead to greater error. The
second method, while more complicated, seems more appropriate and will be
discussed in more detail. It should be noted that this processing is performed
off-line, as real-time operation is not required.

The particular goal of the initial phase of this research project is to
achieve a quantitative estimate of MTs’ interaction when a certain number of
them are presented in a given frame. To clarify the definition of interaction
between MTs, one should understand that two MTs apparently crossing each other
in the time-lapse image does not necessarily mean that they are interacting.
In fact, they might only pass each other without any physical connection at
all. The two discriminating features that help classify between MTs that only
pass each other versus those that interact with each other are their direction
and velocity while crossing. Two MTs just passing each other exhibit no
significant change in their velocity and direction functions, while
interacting MTs will show a completely different behavior. The next step is
feature extraction which requires the exact segmentation of each of the
crossing MTs first.

Segmenting a single MT in a frame that includes several MTs which may be
crossing each other is not a trivial task. To accomplish this, the initial
step is the very common method of thresholding the image (i.e., a single
frame). Having obtained a binary thresholded image, the next step is to apply
the region growing algorithm to find all the neighboring non-zero gray-level
pixels indicating distinct objects. This is accomplished based on a simple
iterative searching algorithm for non-zero values in an 8-pixel wide
neighborhood that starts from an initial pixel. As a result we will get some
detected regions which might belong to a single MT or a crossing situation.
Obviously, in case of a single MT the segmenting mission is already
accomplished, while crossing regions are the target of more processing to
extract single MTs out of a crossing situation. Crossing regions can be
distinguished by their higher than usual gray level (due to the added gray-
level per each MT’s pixels, the crossing area is at least two times“brighter”
than its surrounding pixels). Since this simple sequence of methods
(thresholding and region growing) fails to segment individual MTs out of a
crossing situation, several methods have been explored for this purpose. The
remainder of this subsection describes these methods, along with their cons
and pros.

####  2.2.1 Template Matching

Template-matching is the first applied method to segment single MTs out of
crossed sections. Template matching finds the spatial position of a given
pattern by performing a pixel-wise comparison of the image to the given
template of desired pattern [9]. However, defining an exhaustive collection of
all possible crossing scenarios and picking the most similar template to
properly segment passing MTs seems impractical. The matching criteria can be
defined in various ways: correlation function, squared error, or a simple sum
of pixel by pixel multiplication. Herein, we used 3 lengths of 5, 15, and 50
pixels along with 4 different orientations of:
${\\{{\frac{\pi}{4}},{\frac{\pi}{2}},{\frac{3\pi}{4}},\text{ and }{\pi}\\}}$
to build our pool of templates. Results of applying different criteria
functions on the training data are compared and the best is used to quantify
the test data in the results section.

####  2.2.2 Gradient-based Method

Another method was investigated which is based on the existing contrast in
each object. This method uses the relatively large gradient value at the edges
to produce a binary image indicating boundary pixels. An oriented iterative
searching algorithm then allocates these pixels to distinct objects. This
orientation is specified by gradient vector and helps select between clockwise
or counterclockwise direction of search in the 8-pixel wide neighborhood. The
criteria is having the least gradient change in the extracted results. After a
closed path of boundary pixels, inter-pixels would be determined next.
Although implementing this method is more time intensive, it can work better
in special situations compared to the previously described methods. Herein, we
utilized the training data to measure the threshold value for the gradient
change under which we can discard the changes in the test data.

####  2.2.3 Morphological-based Method

The next method is based on using the morphological erosion of a square-shape
of size $w$ (predetermined width of MTs) on the binary thresholded version of
the frames. Doing so leads to binary images consist of lines or unit-width MTs
called skeletons. Removing the connector pixels (the higher gray-level ones at
cross section) results in some lego pieces. Different combinations of these
pieces and connectors can montage different potential segmentation solutions.
Applying the correlation-based matching criteria (which is the ideal indicator
based on our training data) to the last frame’s results helps to resolve the
uncertainty for segmentation of the object in the current frame. Although the
proposed heuristic algorithm seems to work fine, its practical implementation
fails to improve the results significantly. On the other hand, this method is
not a fully automatic method since it necessitates having the ground truth for
individual MT segmentation in at least one initial frame.

####  2.2.4 Adaptive Template Matching

Another solution is to take the time resolution into account so that the
information in the past can be fused into our segmentation results at the
current frame. This method is similar to the classic template matching.
However, transiting through the frames, we update the pool of templates by
incorporating the segmentation results from the former frame. Since
correlation-based matching criteria presented superior results while
processing the images from the training dataset, we utilized it for segmenting
individual MTs in the test data. This method also is semi-automatic as it
depends on the manual results of individual MT segmentation in at least one
initial frame.

###  2.3 Processing: Data Association in a Tensor-based framework

General implementation of this idea creates a tensor in which, $N$ horizontal
$m\times m$ planes present $N$ frames containing $m$ objects. The $i^{th}$ row
and $j^{th}$ column element at each plane keeps the value of minimum distance
existing between $i^{th}$ and $j^{th}$ MTs of the relevant frame and they can
be assigned to either zero or one depending on the history of their behavior.
Transition to a new frame adds another horizontal plane to the tensor. These
horizontal planes are symmetric and have unit value on their main diagonal.
Sudden appearance and disappearance of any new MT is exposed to the tensor in
form of inserting new vertical planes or zeroing out the relevant row and
column for the next frames. This tensor keeps a history of MTs’ relative
positions toward each other along time. Any of the previously presented
segmentation method followed by the template matching are used to populate the
tensor. Then the same process is applied to all the frames in reverse order so
that the information in upcoming frames can be used to improve the results.
Any dissimilarity in reverse processing indicates further investigation may be
needed. Minimizing the costs of multiple MT association along the sequential
frames lead us to have the optimal extracted trajectories.

###  2.4 Problem Restatement

As it was mentioned in details before, the problem to be solved in this
project is detecting the MT-MT interaction in a series of time-lapsed images.
This description implies a classification problem in which one attempts to
categorize the located crossing-over MTs into two primary groups of MTs just
passing over each other without any interaction and MTs which are physically
interacted to each other through MAPs.

##### Features of Interest:

The discriminating features that help classify between the aforementioned
categories are MTs’ direction and velocity while crossing. The reason for such
feature selection is rooted in the fact that two MTs just passing each other
have a continuous trend in their velocity and direction functions, while
interacting MTs will show a completely different behavior.

##### Feature Extraction:

The next step is feature extraction that includes deriving velocity and
direction of the MTs while crossing each other. This requires the exact
segmentation of each crossing MT.

##### Tracking and Feature Extraction:

After a comprehensive segmentation, a sequence of the sets of pixels for each
MT is obtained. Every set represents a MT at a frame and includes the pixels’
spatial location and their gray-level values at that frame. By following the
frames in order, one can connect a MT’s relevant sets to track its path. Based
on this information, the head block can be identified next. In the case of a
single MT, features are defined in terms of the average moving direction and
speed of the head block’s pixels. In crossing situations, the same features
are extracted with respect to the overlapped area. Assuming a certain pixel,
moving from $(x_{i},y_{i})$ to $(x_{e},y_{e})$ in a transition over two
consecutive frames, its displacement vector is calculated as in Equation 1.

$\overrightarrow{d}=\begin{bmatrix}(x_{e}-x_{i})\\\
(y_{e}-y_{i})\end{bmatrix}$ (1)

Thus, this pixel’s direction and the speed of the movement are calculated
regarding Equation 2 expressions in which ${t_{0}}$ is the frame to frame time
step. Regarding the definition for direction, recall that the reference $x$
and $y$ axis are respectively in downward vertical and left to right
horizontal directions and the origin is located at the very upper-left corner
of the image.

$\theta=\arctan\frac{(x_{e}-x_{i})}{(y_{e}-y_{i})},\quad
V=\frac{\sqrt{(x_{e}-x_{i})^{2}+(y_{e}-y_{i})^{2}}}{t_{0}}$ (2)

##  3 Results

For segmentation, comparisons are made among all the five methods of
thresholding and region growing only, template matching, gradient-based,
erosion-based and the adaptive template matching, each of which followed by
the tensor-based data association method. The simulated time-lapse image are
provided as their input and results in terms of the pixel coordinates are
obtained accordingly. The real data contains 23 RGB, time-lapse sequences with
each having a duration of $23.16\pm 12.68$ seconds. Each sequence is sampled
at 16 frames per second rate with $256\times 256$ pixels spatial dimension. In
order to better present the results, the percentage of the correct segmented
pixels in every frame and for every object has been acquired for each method.
These results are averaged over $20$ frames while tracking $3$ specific MTs
and are depicted in Figure 2. According to each object’s movement pattern and
interactions, different results are obtained. Object $1$ is a single MT that
hardly passes over other MTs during the target frames. The second object
passes along another MT’s length at the latter frames which brings down the
segmentation results especially in case of thresholding and the gradient-based
methods.

Figure 2: Results of segmentation in form of the average percentage of
correctly segmented pixels over 20 frames. Figure 3: Results of feature
extraction in form of the extracted velocity of a single MT and the average of
30 MTs.

The last MT passes through another MT at the mid-frames and detached at the
end. Results for this scenario look compelling. The first two methods
completely fail to segment the object even after passing the crossed section
in contrast to the adaptive template matching method which keeps the history
of the interactions, uses the temporal information and results in the best
among all.

Figure 3 shows the plots of the extracted velocity using all these five
methods comparing to its true value. The upper plot shows the velocity of a
single MT during 50 sequential frames while the bottom one is the average
velocity of all $30$ MTs along the same frames. As it can be seen in the upper
plot, a change in trend of the velocity caused by the simulated interaction is
happening during frame numbers $80$ to $90$ and in this special case, the only
method that segmented the MT with a true estimation of velocity is the tensor-
based algorithm. To have a general perspective, the average plot of the true
velocities of all MTs comes with the average of velocity estimation results
using the mentioned algorithms. Once again, it can be implied that the first
two methods have the largest deviations from the true value of the average
velocity while adaptive template matching method has the least error. However,
average velocity is not a good metric for evaluation since it disguises the
poor performance of an algorithm by averaging, which is not desired especially
since only interacting scenarios are our actual targets.

##  4 Discussion

Although the results seem promising, there are more complicated scenarios
which need to be investigated. A generalized segmentation method should
capture all possible scenarios, such as situations in which a MT has
interactions with more than one MT at different locations, or more than one
interactions occurs at the same place, or even situations where two MTs have a
long overlap along their length. The processing algorithm should consider any
level of interaction at any location and discriminate between being close
versus overlapping. Using temporal resolution information by keeping a history
of the interactions of skeleton images and applying the correlation-based
template matching showed acceptable initial results, but more accurate methods
are needed for more complicated scenarios. The need for improved modeling is
indicated. Regarding the inherent stochastic nature of the MTs’ motility,
applying probabilistic models and estimating methods may yield more reliable
results.

##  5 Conclusions

In order to detect microtubule-microtubule interactions several methods have
been applied in this paper. Although these methods have simple implementations
and provide reasonable results, their deterministic framework shows a
restricted performance for a stochastic problem. Thus, applying probabilistic
models may yield more reliable estimations. Significant progress has been made
to simulate and extract metrics of video microscopy of microtubules. As this
research progresses, it is hoped that it will provide a better understanding
of the underlying mechanics of cancer.

## References

  * [1] J. Challoner, The Cell: A Visual Tour of the Building Block of Life.  University of Chicago Press,2015
  * [2] K. T. Applegate, Quantitative image analysis algorithms for the measurement of cytoskeleton dynamics.  The Scripps Research Institute, 2010.
  * [3] National University of Singapore, Cellular Organization MB Info. https://www.mechanobio.info/topics/cellular-organization/.  Accessed: 2016-09-30.
  * [4] C. Gell, V. Bormuth, G. J. Brouhard, D. N. Cohen, S. Diez, C. T. Friel, J. Helenius, B. Nitzsche, H. Petzold, J. Ribbe, et al., “Microtubule dynamics reconstituted in vitro and imaged by single-molecule fluorescence microscopy,” Methods in cell biology, vol. 95, pp. 221–245, 2010.
  * [5] N.J. Ganem, S.A. Godinho, D. Pellman, “A mechanism linking extra centrosomes to chromosomal instability,” Nature, vol. 460, no. 7252, pp. 278, 2009.
  * [6] J. Avila, Microtubule proteins.  Boca Raton, Fla. : CRC Press, 1990.
  * [7] S. L. Rogers and J. M. Scholey, “Motility assays for microtubule motor proteins,” eLS, 2004.
  * [8] F. Pampaloni and E.-L. Florin, “Microtubule architecture: inspiration for novel carbon nanotube-based biomimetic materials,” Trends in biotechnology, vol. 26, no. 6, pp. 302–310, 2008.
  * [9] C. E. Walczak and S. L. Shaw, “A map for bundling microtubules,” Cell, vol. 142, no. 3, pp. 364–367, 2010.
  * [10] L. Bruno, M. Salierno, D. E. Wetzler, M. A. Despósito, and V. Levi, “Mechanical properties of organelles driven by microtubule-dependent molecular motors in living cells,” PLoS One, vol. 6, no. 4, p. e18332, 2011\.
  * [11] A. Castillo and M. J. Justice, “The kinesin related motor protein, eg5, is essential for maintenance of pre-implantation embryogenesis,” Biochemical and biophysical research communications, vol. 357, no. 3, pp. 694–699, 2007.
  * [12] Y. Chen and W. O. Hancock, “Kinesin-5 is a microtubule polymerase,” Nature communications, vol. 6, 2015.
  * [13] I. Smal, K. Draegestein, N. Galjart, W. J. Niessen, and E. Meijering, “Particle filtering for multiple object tracking in dynamic fluorescence microscopy images: Application to microtubule growth analysis,” IEEE Transactions on Image Processing, vol. 27, no. 6, pp. 789–804, 2008.
  * [14] M.R. Dilsaver, M. Tomschik, C. Geisler,J.C. , “Tau-based fluorescent protein fusions to visualize microtubules,” Cytoskeleton, vol. 74, no. 6, pp. 221–232, 2017.


# Effectiveness
of the Bendixson–Dulac theorem

Armengol Gasull  Departament de Matemàtiques, Edifici Cc, Universitat Autònoma
de Barcelona, 08193 Cerdanyola del Vallès (Barcelona), Spain.  Centre de
Recerca Matemàtica, Edifici Cc, Campus de Bellaterra, 08193 Cerdanyola del
Vallès (Barcelona), Spain<EMAIL_ADDRESS>and  Hector Giacomini  Institut
Denis Poisson. Université de Tours, C.N.R.S. UMR 7013. 37200 Tours, France.
<EMAIL_ADDRESS>
###### Abstract.

We illustrate with several new applications the power and elegance of the
Bendixson–Dulac theorem to obtain upper bounds of the number of limit cycles
for several families of planar vector fields. In some cases we propose to use
a function related with the curvature of the orbits of the vector field as a
Dulac function. We get some general results for Liénard type equations and for
rigid planar systems. We also present a remarkable phenomenon: for each
integer $m\geq 2,$ we provide a simple $1$-parametric differential system for
which we prove that it has limit cycles only for the values of the parameter
in a subset of an interval of length smaller that $3\sqrt{2}(3/m)^{m/2}$ that
decreases exponentially when $m$ grows. One of the strengths of the results
presented in this work is that although they are obtained with simple
calculations, that can be easily checked by hand, they improve and extend
previous studies. Another one is that, for certain systems, it is possible to
reduce the question of the number of limit cycles to the study of the shape of
a planar curve and the sign of an associated function in one or two variables.

###### Key words and phrases:

Limit cycle; Periodic orbit; Bendixson–Dulac theorem; Liénard equation

###### 2010 Mathematics Subject Classification:

Primary: 34C07. Secondary: 37C27

##  1\. Introduction

Despite all the efforts dedicated to solve the second part of the Hilbert’s
16th problem, it is yet a very difficult task to obtain criteria that give
explicit upper bounds for many concrete families of planar smooth vector
fields. Although there is no any universal approach, the aim of this paper is
to present several families of planar systems for which the Bendixson–Dulac
theorem allows to get, in a fast and elegant way, an upper bound of their
number of limit cycles. We will avoid results based on cumbersome
computations.

The families that we will consider include extensions of Liénard systems and
rigid systems. As we will see, we obtain new results and we also present
simple proofs of some recent results. They give explicit upper bounds for
several families of planar vector fields. These upper bounds are also
sharpened when we deal with more particular systems, obtaining results of at
most two, one, or none limit cycles.

Our main results for Liénard type systems are contained in Section 3. They are
given in Theorem 3.1, that deals with a version of Wilson Liénard systems
which always have an algebraic limit cycle, in Theorem 3.4 that studies a
family recently introduced in [26], in Theorem 3.9 that extends the classical
theorem of Massera, and in Theorem 3.7. In fact, this last result includes the
remarkable phenomenon highlighted in the abstract: the family

$\begin{cases}\dot{x}=y-\lambda|y|^{m}(x^{3}-x),\\\ \dot{y}=-x,\end{cases}$

introduced in [26], has for $m\geq 2$, limit cycles only for some values of
$\lambda$ contained in the interval of length $3\sqrt{2}(3/m)^{m/2},$ centered
at the origin. Notice that for $m$ big it is extremely thin. This interval
decreases exponentially with $m$.
Our main result for general rigid systems is given in Theorem 4.1 of Section
4. It is applied to recover, in a simple way, known results for rigid cubic
systems and to a family containing non polynomial vector fields.

It seems to us that not all the mathematical community that works on these
topics is aware of the capability of the Bendixson–Dulac approach. The goal of
this work is double: we try to change this perception and we also present
several new results and easy proofs of some known results. For instance, in
most textbooks, the proof of the uniqueness and hyperbolicity of the limit
cycle for the classical van der Pol equation needs some work. By using this
approach there are extremely simple proofs, see Corollary 3.5 and Remark 3.6.

The today known as Bendixson–Dulac theorem was first formulated by Ivar
Bendixson in 1901 ([1]), and later developed by Henri Dulac in 1937 ([12]). He
improved Bendixson’s result by introducing a new parameterization of the time,
via the today called Dulac functions. This result appears, under different
versions, in most differential equations textbooks. One of the pioneers to try
to go further with this approach was Yamato ([29]). Afterward, one of its main
defenders was Cherkas, who used and developed it, see for instance [4, 5]. The
authors of this work also often apply and try to extend this method, see [15,
16, 17]. More examples about its applicability can be seen in the survey [18].

In this paper we will use the version of the Bendixson–Dulac theorem that we
state below, after introducing some notations and definitions. For
completeness, in Section 2 we present a proof based on the version of
Bendixson–Dulac theorem for multiply connected regions that is proved for
instance in [4, 15, 23].

Given an open connected subset $\mathcal{U}\subset\mathbb{R}^{2},$ with
finitely many holes, we will denote by $\ell=\ell(\mathcal{U})$ this number of
holes, that is, the number of bounded components of
$\mathbb{R}^{2}\setminus\mathcal{U}.$ Notice that if $\mathcal{U}$ is simply
connected then $\ell(\mathcal{U})=0.$ We also set
$\mathbb{R}^{+}=\\{x\in\mathbb{R}\,:\,x>0\\}.$

For a continuous function $f:\mathbb{R}^{n}\to\mathbb{R},$ not changing sign
and vanishing on a null measure set, we will denote by
$\operatorname{sign}(f)$ the sign of $f$ at any of its point where it is not
zero. Moreover, given an equilibrium point or a periodic orbit, when we say
that its stability is given by the sign of $f$ we mean that the object is an
attractor (resp. a repeller) whenever $\operatorname{sign}(f)<0$ (resp.
$\operatorname{sign}(f)>0$).

######  Definition 1.1.

Given a function $V:\mathbb{R}^{2}\to\mathbb{R}$ of class $\mathcal{C}^{1}$ we
will say that it is admissible if:

  1. (i)

The vector $\nabla V$ vanishes on $\\{V(x,y)=0\\}$ at finitely many points.

  2. (ii)

The set $\\{V(x,y)=0\\}$ has finitely many connected components.

  3. (iii)

The set $\mathbb{R}^{2}\setminus\\{V(x,y)=0\\}$ has $j$ connected components,
$\mathcal{U}_{i},i=1,2,\ldots j,$ and for all of them
$\ell(\mathcal{U}_{i})<\infty.$

Associated to $V,$ we define the non negative integer number

$L(V):=\sum_{i=1}^{j}\ell(\mathcal{U}_{i}).$

######  Theorem 1.2 (Bendixson–Dulac theorem).

Consider a $\mathcal{C}^{1}$ planar differential system

$\dot{x}=P(x,y),\quad\dot{y}=Q(x,y),$ (1)

and denote by $X=(P,Q)$ its associated vector field. Let
$V:\mathbb{R}^{2}\to\mathbb{R}$ be an admissible function such that there
exists $s\in\mathbb{R}^{+}$ for which the function

$M_{s}:=\frac{\partial V}{\partial x}P+\frac{\partial V}{\partial
y}Q-s\left(\frac{\partial P}{\partial x}+\frac{\partial Q}{\partial
y}\right)V$ (2)

does not change sign and vanishes only on a null measure set. Define

$L_{X}(V):=N+L(V),$

where $N$ is the number of periodic orbits of $X$ contained in the set
$\mathcal{V}=\\{V(x,y)=0\\}.$

Then, the differential system (1) has at most $L_{X}(V)$ periodic orbits,
which are limit cycles. Moreover, each limit cycle not contained in
$\mathcal{V}$ is hyperbolic, it is contained in one of the connected
components $\mathcal{U}_{i}$ of $\mathbb{R}^{2}\setminus\mathcal{V}$ and, for
each $i=1,2,\ldots,j,$ there are at most $\ell(\mathcal{U}_{i})$ limit cycles
in the component $\mathcal{U}_{i}.$ The stability of each of these limit
cycles is given by the sign of $-VM_{s}$ on the region $\mathcal{U}_{i}$.

######  Remark 1.3.

The function $M_{s},$ when $s\leq 0,$ can also be used to control te number of
limit cycles of (1), see [4, 18]. In particular, notice that $M_{0}=\dot{V}$
and it can be readily seen that, when $s<0,$ the theorem also works, giving
that $L_{X}(V)=N.$ In this work, we do not use this range of values of $s.$ In
fact, in most of our applications we will use $s=1,$ although the values $s=2$
and $s=1/3,$ also will appear.

Observe also that, somehow, this version of the Bendixson–Dulac theorem
relates the second part of the Hilbert’s 16th problem, which deals with the
number of limit cycles ([22]), with the first part, that deals with the number
and distribution of the ovals of a planar algebraic curve ([27]).

Similarly of what happens when one tries to use Lyapunov functions, the main
difficulty in the above theorem for its practical use is the choice of the
function $V$ and of the positive real number $s.$ In other words, the choice
of a suitable Dulac function. As we will see, the function that gives the
curvature of the orbits of (1) is sometimes a good candidate for $V.$

Moreover, the most difficult condition to be checked is that $M_{s}$ does not
change sign. Hence, several approaches try to arrive to functions for which
this question can be more easily studied. For instance, one of these
situations is when it is a function of only one variable or the product of two
functions of one variable, see again [18] for some examples. Another one is
when, from some point of view, we can look to $M_{s}$ as a quadratic
polynomial.

Finally, notice that, given $V$ and $X,$ the computation of the number $N$ in
$L_{X}(V)$ is usually not difficult, while $L(V)$ depends on the topology of
the set $\\{V(x,y)=0\\},$ see Section 2.2. When $V$ is polynomial in one of
its variables, to get an upper bound of $L(V)$ is an affordable task.

##  2\. Preliminary results

For the sake of notation, from now on, in this paper we will denote the
partial derivatives as subscripts. Hence, for instance, for $F=F(x,y),$
$F_{x}=\frac{\partial F}{\partial x},$ or
$F_{x,y}=\frac{\partial^{2}F}{\partial x\partial y}$.

###  2.1. Proof of Theorem 1.2

We recall a version of the Bendixson–Dulac theorem for multiply connected
regions, see [4, 15, 23].

######  Theorem 2.1.

Consider a $\mathcal{C}^{1}$ planar differential system

$\dot{x}=P(x,y),\quad\dot{y}=Q(x,y),$ (3)

defined on $\mathcal{U}\subset\mathbb{R}^{2},$ an open subset such that
$\mathbb{R}^{2}\setminus\mathcal{U}$ has $\ell$ bounded components, and denote
by $X=(P,Q)$ its associated vector field. Let
$B:\mathcal{U}\to\mathbb{R}^{+}\cup\\{0\\}$ be a $\mathcal{C}^{1}$ function
such that

$\operatorname{div}(BX)=(BP)_{x}+(BQ)_{y}$

does not change sign and vanishes only on a null measure set. Then the system
(3) has at most $\ell$ limit cycles in $\mathcal{U}.$ Moreover, all of them
are hyperbolic and their stability is given by the sign of
$\operatorname{div}(BX).$

To prove Theorem 1.2 first one has to show that the periodic orbits of (1) are
either contained in $\mathcal{V}$ or do not cut this set. This fact follows
because $M_{s}\big{|}_{\mathcal{V}}=\nabla V\cdot X=\dot{V}$ does not change
sign. Hence system (1) can have some periodic orbits contained in
$\mathcal{V},$ say that it has $N,$ and all the others that are strictly
contained in each of the connected components of
$\mathbb{R}^{2}\setminus\mathcal{V}.$ Fix one of these connected components,
say $\mathcal{U}_{i}.$ To control the number of periodic orbits in this set we
will apply Theorem 2.1 with $B=|V|^{-1/s}$ and $\mathcal{U}=\mathcal{U}_{i}.$
Some computations give that

$\operatorname{div}\Big{(}|V|^{-1/s}X\Big{)}=-\frac{1}{s}\operatorname{sign}(V)|V|^{-1/s-1}M_{s}$
(4)

and, by hypothesis, this function does not change sign on $\mathcal{U}_{i}.$
As a consequence, the maximum number of periodic orbits in $\mathcal{U}_{i}$
is $\ell(\mathcal{U}_{i}),$ as we wanted to prove. Moreover, by using (4) and
again Theorem 2.1 we get that all of them are hyperbolic and their stability
is given by the sign of $-VM_{s}.$

######  Remark 2.2.

Notice that in Theorem 1.2 nothing is said about the hyperbolicity of the
limit cycles contained in $\mathcal{V}.$ As we will see in Corollary 3.2, they
can be hyperbolic or not.

###  2.2. About the practical calculation of $L(V)$

Given an admissible function $V,$ the computation of $L(V)$ relies on the
study of the topology of each of the connected components $\mathcal{U}_{i},$
of $\mathbb{R}^{2}\setminus\mathcal{V},$ where $\mathcal{V}=\\{V(x,y)=0\\}.$
Then $L(V)$ is the sum of all the quantities $\ell(\mathcal{U}_{i})$, where
these values are the number of bounded components of
$\mathbb{R}^{2}\setminus\mathcal{U}_{i}.$ In fact, it also holds that the
fundamental group of $\mathcal{U}_{i}$ is
$\pi_{1}(\mathcal{U}_{i})=\mathbb{Z}\stackrel{{\scriptstyle\ell)}}{{*\cdots*}}\mathbb{Z},$
where $\ell=\ell(\mathcal{U}_{i}).$ In all concrete situations appearing in
this work there is a more direct way for obtaining $L(V).$ This number is
simply the number of bounded connected components of $\mathcal{V}.$

###  2.3. Curvature of the orbits

It is know that the function

$K^{\perp}:=Q^{2}P_{x}+P^{2}Q_{y}-PQ(P_{y}+Q_{x}),$

that is the numerator of the curvature of the orbits of the vector field
$X^{\perp}=(-Q,P),$ orthogonal to the vector field $X=(P,Q),$ associated to
the system (1), can be used to know the stability of the periodic orbits of
(1) and other dynamical features of its phase portrait, see [6, 10, 14] or
[30, p. 29]. For instance, Diliberto in 1950 proved that a limit cycle is
hyperbolic and stable (resp. unstable) if and only if

$\int_{0}^{l}K^{\perp}(\gamma(s))\,{\rm d}s<0\quad(\mbox{resp.}\,\,>0),$

where $\gamma(s)$ is its parameterization by the arc length and $l$ is its
length.

In this work we will see that the function

$K:=Q^{2}P_{y}-P^{2}Q_{x}+PQ(P_{x}-Q_{y}),$ (5)

proportional to the numerator of the curvature of the orbits of $X$ is, in
several cases, a good candidate for a suitable choice of $V$ in Theorem 1.2.
Notice that $K=Q\dot{P}-P\dot{Q}=Q(P_{x}P+P_{y}Q)-P(Q_{x}P+Q_{y}Q).$ As far as
we know, this is the first time that this function $K$ is used to control the
number of limit cycles of planar differential systems. We prove:

######  Theorem 2.3.

Consider planar system (1) of class $\mathcal{C}^{2}.$ Assume that the
function

$D:=P^{2}Q\big{(}P_{xx}-2Q_{xy}\big{)}+PQ^{2}\big{(}2P_{xy}-Q_{yy}\big{)}+Q^{3}P_{yy}-P^{3}Q_{xx}$

does not change sign and vanishes on a null measure set. Then the system (1)
has at most $L_{X}(V)$ limit cycles, where $V=K$ is given in (5) and
$L_{X}(V)$ is defined in Theorem 1.2.

###### Proof.

By taking $V=K,$ as in equation (5), and $s=1,$ the function $M_{1}$ given in
Theorem 1.2 is $M_{1}=D$ and the theorem follows. ∎

We will apply this result at the end of Section $3$ for Liénard systems and in
Section $4$ to rigid systems.

##  3\. Liénard type systems

We present several applications of the Bendixson–Dulac theorem to two families
related with Liénard systems.

###  3.1. Liénard systems with an explicit solution

We study a family of Liénard type equations introduced recently in [21] that
includes the Wilson family of Liénard equations ([28]), which gave the first
example of such equations having an algebraic limit cycle. More concretely, we
consider systems

$\begin{cases}\dot{x}=y-(x^{2}-1)B(x),\\\ \dot{y}=-x(1+yB(x)),\end{cases}$ (6)

where $B$ is a $\mathcal{C}^{1}$ function. They have the invariant algebraic
curve $C(x,y)=x^{2}+y^{2}-1=0,$ because $C_{x}P+C_{y}Q=-2xBC,$ where $X=(P,Q)$
denotes the vector field associated to (6). Hence, when this system has not
equilibrium points on the curve, it is a periodic orbit. Moreover, depending
on the choice of the function $B$, it can be a limit cycle.

The following result allows to extend, and to prove in an easier way, the
recent results about the maximum number of limit cycles of the above system
when $B(x)=x^{3}-bx$ given in [3, 21].

######  Theorem 3.1.

Consider the system (6) with $B(x)=x\int_{0}^{x}W(t)/t\,{\rm d}t-bx,$ where
$W$ is any function that does not change sign, vanishes at isolated points,
and such that $B$ is of class $\mathcal{C}^{1}.$ Then this system has at most
$L+N$ limit cycles, where $L$ is the number of bounded connected components of
the set $\mathcal{B}=\\{x\in\mathbb{R}\,:\,(B(x)+2x)(B(x)-2x)\geq 0\\}$ plus
one, and $N\in\\{0,1\\}.$ In fact $N=1$ when
$\mathcal{C}=\\{x^{2}+y^{2}-1=0\\}$ is free of equilibrium points of the
system, and then this set is one of the limit cycles, and $N=0$ otherwise.
Moreover, all the limit cycles but $\mathcal{C}$ are hyperbolic and their
stability is given by the sign of $VW$ in the connected component of
$\mathbb{R}^{2}\setminus\\{V(x,y)=0\\}$ where they lie, with

$V=(1-x^{2}-y^{2})\left(x^{2}+y^{2}+B(x)y\right).$ (7)

###### Proof.

Consider the function $V$ given in (7) and $s=1$ in Theorem 1.2. Then,

$M_{1}(x,y)=x(x^{2}+y^{2}-1)^{2}\big{(}B(x)-xB^{\prime}(x)\big{)}=-x^{2}(x^{2}+y^{2}-1)^{2}W(x).$

Hence, thanks to the imposed conditions on $W,$ we can apply Theorem 1.2.
Moreover, since $\mathcal{C}$ is invariant, and contained in
$\mathcal{V}=\\{V(x,y)=0\\},$ we have that $N\in\\{0,1\\}$ and the number of
limit cycles of the system is bounded by $L(V)+N.$ To get $L(V)$ we study the
bounded connected components of $\mathcal{V},$ see Section 2.2. Notice that
these components are formed by the oval $\mathcal{C}$ together with the
bounded connected components of $x^{2}+y^{2}+B(x)y=0.$ Since this curve also
writes as

$y=\frac{-B(x)\pm\sqrt{(B(x)+2x)(B(x)-2x)}}{2},$

it is clear that these components are obtained by joining the curves plus and
minus defined for $x$ on each of the bounded connected components of
$\mathcal{B}.$ Hence $L(V)$ is at most $L$ and the theorem follows. ∎

The following corollary gives an easier and different proof of all the results
about the maximum number of limit cycles of (6) when $B(x)=x^{3}-bx,$

$\begin{cases}\dot{x}=y-(x^{2}-1)(x^{3}-bx),\\\
\dot{y}=-x\big{(}1+y(x^{3}-bx)\big{)},\end{cases}$ (8)

obtained in [3, 21]. It also solves in the best possible way the some times
called Coppel’s problem for polynomial systems, which in his own words (when
restricted to quadratic systems) says:“Ideally one might hope to characterize
the phase portraits of quadratic systems by means of algebraic inequalities on
the coefficients,” see [9]. The relevant values describing the bifurcations of
the limit cycles of this system are $\underline{b}$ and $b^{*},$ see again
[21]. We will show below that the value $\underline{b}\approx-1.44$ is
algebraic. It is the negative root of the polynomial

$4b^{6}-12b^{5}-4b^{4}+28b^{3}+56b^{2}-72b-229=0,$ (9)

which is invariant under the change of variables $b\rightarrow 1-b.$ The
quantity $b(1-b)$ satisfies a third degree polynomial equation and then it is
possible to express all the roots in terms of radicals but we prefer to omit
their explicit expressions because they are rather complicated. The value
$b^{*}\approx 0.747$ is the only zero of the function
$Z:(\underline{b},1-\underline{b})\to\mathbb{R},$

$Z(b)=\int_{0}^{1}\frac{8(b-3x^{2})\sqrt{1-x^{2}}}{x^{8}-(2b+1)x^{6}+(b+2)bx^{4}-b^{2}x^{2}+1}\,{\rm
d}x.$

The sign of this function gives the stability of
$\mathcal{C}=\\{x^{2}+y^{2}-1=0\\}$ when $\mathcal{C}$ is a limit cycle. Most
probably $b^{*}$ is a non-algebraic number. The function $Z$ was obtained in
[21] from the integral of the divergence of the system on the algebraic limit
cycle after some algebraic manipulations. In fact, the discriminant with
respect to $x$ of the denominator of the integrand gives the polynomial of the
left-hand side of (9) that determines $\underline{b}.$

######  Corollary 3.2.

System (8) has at most two limit cycles, taking into account their
multiplicities. More concretely:

  1. (i)

It has no limit cycle for $b\leq\underline{b}.$

  2. (ii)

$\mathcal{C}$ is its only limit cycle, and it is hyperbolic and attractor when
$b\in(\underline{b},0].$

  3. (iii)

It has two limit cycles, one hyperbolic, repeller and surrounded by
$\mathcal{C},$ and $\mathcal{C}$ itself, which is hyperbolic and attractor,
when $b\in(0,b^{*}).$

  4. (iv)

$\mathcal{C}$ is its only limit cycle, and it is double and semi-stable when
$b=b^{*}.$

  5. (v)

It has two limit cycles, one hyperbolic, attractor and surrounding
$\mathcal{C},$ and $\mathcal{C}$ itself, which is hyperbolic and repeller,
when $b\in(b^{*},1-\underline{b}).$

  6. (vi)

It has one limit cycle surrounding $\mathcal{C}$ that is hyperbolic and
attractor, when $b\geq 1-\underline{b}.$

###### Proof.

First, we will prove the most difficult part: the maximum number of limit
cycles of the system is three. This will essentially be a direct consequence
of Theorem 3.1. All the other results about this system can be obtained from
the standard techniques of the qualitative theory of planar differential
systems.

When $b\leq 0$ the only limit cycle is $\mathcal{C}$ because in polar
coordinates, $\dot{r}=r(r^{2}-1)(b-r^{2}\cos^{2}\theta)\cos^{2}\theta$ does
not vanish outside $\mathcal{C}.$

When $b\geq 3/2,$ the maximum number of limit cycles is two. To prove this we
apply Theorem 1.2 with $V(x,y)=x^{2}+y^{2}-1$ and $s=1/3.$ Then

$M_{1/3}(x,y)=\frac{1}{3}(x^{2}+y^{2}-1)\big{(}(2b-3)x^{2}+b\big{)}\geq 0.$

Since for these values of $b,$ $M_{1/3}$ does not vanish outside $\mathcal{C}$
the maximum number of limit cycles is two, one being $\mathcal{C}$ and at most
another one can exist, and in this case it must surround $\mathcal{C}.$

Finally, consider the values of $b\in(0,3/2).$ In fact, we can consider
$b\in[0,2].$ We apply Theorem 3.1 with $W(x)=2x^{2}\geq 0.$ We get that
$B(x)=x^{3}-bx,$ and

$\displaystyle\mathcal{B}$
$\displaystyle=\\{x\in\mathbb{R}\,:\,x^{2}(x^{2}-2-b)(x^{2}+2-b)\geq 0\\}$
$\displaystyle=\big{(}-\infty,-\sqrt{2+b\,}\,\big{]}\cup\\{0\\}\cup\big{[}\sqrt{2+b\,},\infty\big{)}.$

Hence, the number of bounded connected components of $\mathcal{B}$ is one and,
as a consequence, $L(V)=2$ and the maximum number of limit cycles is three.
Also, from the proof we know that if the three limit cycles exist, one is
$\mathcal{C},$ there is at most another one, say $\gamma,$ surrounded by
$\mathcal{C},$ and a third one $\Gamma,$ surrounding $\mathcal{C}$.

To reduce this upper bound of three limit cycles by one it suffices to
consider the stability of the origin, the infinity, the possible limit cycles
and the invariant set $\mathcal{C}.$ In fact we have that,

  1. (A)

The stability of the origin is given by the sign of $b.$ Moreover, it is not
difficult to see that when $b\in(\underline{b},1-\underline{b})$ the origin is
the only equilibrium point of the system and that, otherwise, there are also
other equilibrium points, but all of them are on $\mathcal{C}.$

  2. (B)

The set $\mathcal{C},$ which is always invariant by the flow, is a limit cycle
if and only if $b\in(\underline{b},\overline{b}).$ Moreover it is hyperbolic
and stable if $b\in(\underline{b},b^{*}),$ hyperbolic and repeller if
$b\in(b^{*},\overline{b}),$ and semi-stable and double when $b=b^{*}.$ In
fact, in this later case it is repeller from its interior and attractor from
its exterior, see [21]. Moreover, it is also proved in that paper, that when
$b\geq\overline{b},$ the set $\mathcal{C},$ that it is no more a periodic
orbit, is also a repeller.

  3. (C)

The infinity is repeller for $b>0,$ see again [21].

  4. (D)

For $b\in(0,2),$ whenever they exist, $\gamma$ is hyperbolic and repeller and
$\Gamma$ is hyperbolic and attractor. This is a consequence of Theorem 3.1,
because their respective stabilities are controlled by the sign of $VW,$ that
coincides with the sign of $1-x^{2}-y^{2},$ because

$VW=2x^{2}(1-x^{2}-y^{2})(x^{2}+y^{2}+(x^{3}-bx)y),$

and for these values of $b$ the limit cycles must lie in the region
$\\{x^{2}+y^{2}+(x^{3}-bx)y>0\\}$ because it is the only connected component
$\mathcal{U}$ of $\mathbb{R}^{2}\setminus\mathcal{V}$ with
$\ell(\mathcal{U})\neq 0.$

  5. (E)

For $b\geq 3/2,$ $\gamma$ never exits and $\Gamma$ is also hyperbolic and
atractor, because as we have proved above by using Theorem 1.2, its stability
is also given by the sign of

$\qquad-V(x,y)M_{1/3}(x,y)=-\frac{1}{3}(x^{2}+y^{2}-1)^{2}\big{(}(2b-3)x^{2}+b\big{)}\leq
0.$

For instance we will prove item $(v).$ All the other cases follow similarly.
First notice that by (B), $\mathcal{C}$ is a hyperbolic and repeller limit
cycle. Recall that we already have proved that the system has at most one
limit cycle surrounded by $\mathcal{C},$ and another one surrounding
$\mathcal{C}.$ Moreover, whenever they exist they are hyperbolic and their
stabilities are given in (D). By (A) and (C), since the origin is attractor
and the infinity is repeller, we get that there is no limit cycle surrounded
by $\mathcal{C}$ and there is exactly one limit cycle, hyperbolic and stable,
surrounding $\mathcal{C}.$

∎

######  Remark 3.3.

System (8) can be transformed into the classical Liénard system

$\begin{cases}\dot{x}=y-bx+x^{3}+\frac{4b}{3}x^{3}-\frac{6}{5}x^{5},\\\
\dot{y}=-x+b^{2}x^{3}-b(2+b)x^{5}+(1+2b)x^{7}-x^{9}.\end{cases}$ (10)

By using Theorem 1.2 with $s=1$ and $V(x,y)=A(x,y)B(x,y)$ where

$\displaystyle A(x,y)=$
$\displaystyle-225+225x^{2}+25b^{2}x^{6}-30bx^{8}+9x^{10}$
$\displaystyle+(150bx^{3}-90x^{5})y+225y^{2},$ $\displaystyle B(x,y)=$
$\displaystyle 225x^{2}-75b^{2}x^{4}+5b(24+5b)x^{6}-15(3+2b)x^{8}+9x^{10}$
$\displaystyle+(-225bx+25(9+6b)x^{3}-90x^{5})y+225y^{3},$

we get that $M_{1}=2x^{4}A^{2}(x,y)\geq 0.$ Hence, in these variables an upper
bound of the number of limit cycles of system (8) can also be obtained. This
example illustrates that although sometimes it is difficult to find a suitable
$V$ to apply Theorem 1.2, it seems to exist.

###  3.2. Some extended Liénard systems

We consider planar differential equations of the form

$\begin{cases}\dot{x}=y-|y|^{m}F(x),\\\ \dot{y}=-G^{\prime}(x)/2,\end{cases}$
(11)

where $F$ and $G^{\prime}$ are $\mathcal{C}^{1}$ functions satisfying $F(0)=0$
and $G(x)=x^{2k}+o(x^{2k}),$ $m\in\mathbb{N}\cup\\{0\\}$ and $k\in\mathbb{N}.$
Notice that when $G(x)=x^{2}$ and $m=0,$ they include the classical second
order Liénard equations $\ddot{x}+F^{\prime}(x)\dot{x}+x=0.$ The factor
$|y|^{m}$ is added following the recent work [26], where this interesting
system was studied for the first time. Notice that if instead of
$y-|y|^{m}F(x)$ we consider the same system but with the first component equal
to $y-y^{m}F(x),$ then, when $m$ is odd, it would be invariant by the change
of variables and time $(x,y,t)\to(x,-y,-t)$ and the origin would be a
reversible center.

In all our study we skip the case $m=1,$ where the associated vector field is
not of class $\mathcal{C}^{1}.$ In any case, for $m=1,$ and on each of the
regions $y>0$ and $y<0,$ the vector field is integrable (it corresponds to a
differential equation of separated variables) and by using the level curves of
the corresponding first integrals, their phase portraits are easier to be
studied. This approach is the one used in [26] for this case, when
$G(x)=x^{2}.$

######  Theorem 3.4.

Consider the differential system (11) with $m\neq 1.$ If the function
$H:=(m-1)FG^{\prime}+2F^{\prime}G$ does not change sign and vanishes at
isolated points, then the system has at most $J$ limit cycles, all of them
hyperbolic, where $J$ is the number of zeroes of $G^{\prime}.$ In particular,
if $G^{\prime}$ only vanishes at the origin the differential system has at
most one limit cycle.

###### Proof.

We apply the Bendixson–Dulac theorem with $V(x,y)=G(x)+y^{2}-y|y|^{m}F(x)$ and
$s=1.$ Simple computations give that

$M_{1}=\frac{1}{2}|y|^{m}H(x).$

Therefore, since $M_{1}$ satisfies the hypothesis of the Bendixson–Dulac
theorem we have already proved that system (11) has at most $L_{X}(V)$ limit
cycles. We claim that $L_{X}(V)\leq J.$ Since the set
$\mathcal{V}=\\{V(x,y)=0\\}$ does not contain solutions of the differential
system the claim will follow if we prove that $\mathcal{V}$ has at most $J$
bounded connected components, see Section 2.2. Notice that each of these
components can be an oval, an isolated point, or a more complicated set.

To prove this last assertion we first count the number of points of
$\mathcal{V}\cap\\{x=x_{0}\\},$ taking into account their multiplicity, and we
call it $K(x_{0}).$ When $m=0$ it is clear that $K(x_{0})\leq 2,$ because
$V(x_{0},y)=0$ is a quadratic equation in $y.$ When $m\geq 2,$ the equation
$V(x_{0},y)=0$ splits into two trimonomial equations, one for $y\geq 0$ and
another one for $y\leq 0.$ By applying the Descarte’s rule of signs to both
equations, since the monomial $y^{2}$ appears in both, it can be seen that
$K(x_{0})\leq 3.$

Notice that since on $\mathcal{V},$ $M_{1}=\dot{V},$ each bounded connected
component of $\mathcal{V}$ delimits some region either positively or
negatively invariant, and as a consequence its interior must contain at least
one equilibrium point $(x^{*},y^{*})$ of the system. Observe also that even
when the system has other equilibrium points on the line $\\{x=x^{*}\\},$ only
one connected component of $\mathcal{V}$ can cut this line, because
$K(x^{*})\leq 3.$ Hence the bounded connected components of the set
$\mathcal{V}$ must cut the lines $\\{x=x^{*}\\},$ where $G^{\prime}(x^{*})=0,$
and at most one of them cuts each of the lines. As a consequence,
$\mathcal{V}$ has at most $J$ bounded connected components, and $L(V)\leq J$
as we wanted to prove. ∎

Theorem 3.4 can be applied to several differential systems (11). For $m$ and
$G$ fixed, let $W$ be a function that does not change sign, vanishes at
isolated points, and such that the initial value problem for the linear
differential equation

$(m-1)F(x)G^{\prime}(x)+2F^{\prime}(x)G(x)=W(x),\quad F(0)=0,$ (12)

has a regular solution $F.$ Notice that (12) is singular at the zeroes of $G$
and we impose that $F$ must be smooth at these points. Then the correspondent
differential system (11) is under the hypotheses of the theorem. By using this
point of view, we have obtained several families of differential systems where
it is easy to impose that their corresponding functions $H$ do not change sign
and, as a consequence, Theorem 3.4 can be applied. We will skip all the
hypotheses that must be added to guarantee the desired property for $H,$ and
the other ones that the functions $F$ and $G$ must satisfy to fulfill all the
other hypotheses of the theorem, because the reader can easily figure out
them. These families are:

  1. (i)

When $F(x)=A^{p}(x)A^{\prime}(x)B(x),$
$G(x)=cA^{q}(x)(A^{\prime}(x))^{2}B^{2}(x),$ and $m=0.$ Then, it holds that

$H(x)=(2p-q)cA^{p+q-1}(x)(A^{\prime}(x))^{4}B^{3}(x).$

  2. (ii)

When $m=0,$

$\displaystyle F(x)$
$\displaystyle=A^{2p}(x)A^{\prime}(x)B^{q+1}(x),\,\,\,\mbox{and}$
$\displaystyle G(x)$ $\displaystyle=cA^{4p}(x)(A^{\prime}(x))^{2}B^{q}(x),$

we get that

$H(x)=(q+2)cA^{6p}(x)B^{2q}(x)(A^{\prime}(x))^{3}B^{\prime}(x)$

  3. (iii)

When $m=2k,$ $G(x)=x^{2k},$ and

$F(x)=\frac{1}{2}x^{k(1-2k)}\int_{0}^{x}y^{k(2k+1)}Z(y)\,{\rm d}y,$

we obtain that $H(x)=x^{4k}Z(x).$

  4. (iv)

When $m=0,$ $F(x)=a(x^{3}/3-x)$ and

$G(x)=x^{2}-\Big{(}\frac{a^{2}}{8}+6b\Big{)}x^{4}+\Big{(}\frac{a^{2}}{48}+b\Big{)}x^{6},$

we get that

$H(x)=\frac{a(16-3a^{2}-144b)}{12}x^{4}.$

Now we will study in more detail some particular sub-cases of the above
families and we will refine the upper bound for their number of limit cycles
given in Theorem 3.4.

We start with an example contained in the family given in item (i). It
corresponds to $p=1,$ $q=0,$ $c=1/4,$ $A(x)=x^{3}/3-x^{2}/2$ and $B(x)=-2,$
and writes as

$\begin{cases}\dot{x}=y+\dfrac{1}{3}x^{3}(x-1)(2x-3),\\\\[5.69046pt]
\dot{y}=-x(x-1)(2x-1),\end{cases}$ (13)

with $G(x)=x^{2}(1-x)^{2}.$ We will prove that this system has at most one
limit cycle, hyperbolic and stable. The existence of this limit cycle, that
surrounds the three equilibrium points of (13), can be numerically confirmed.

By using Theorem 3.4 when $m=0$ and with

$V(x,y)=x^{2}(x-1)^{2}+\frac{1}{3}x^{3}(x-1)(2x-3)y+y^{2}$

we get that $H(x)=-4x^{4}(x-1)^{4}<0$ and the maximum number of limit cycles
of the corresponding system is three, which is the number of zeroes of
$G^{\prime}(x)=2x(x-1)(2x-1).$ This upper bound can be reduced to two studying
in more detail the set $\mathcal{V}.$ This set is formed by two isolated
critical points located at $(0,0)$ and $(1,0)$, and two disjoint curves going
from infinity to infinity. The point $(0,0)$ is a weak focus and the point
$(1,0)$ is a strong stable focus. The third critical point, located at
$(1/2,-1/24)$, is a saddle point. By computing the first Lyapunov quantity
associated to the weak focus at the origin we conclude that this point is
repulsive. In fact, $\mathbb{R}^{2}\setminus\mathcal{V}$ is formed by three
open sets, two are simply connected and the third one has two holes (the two
critical points located on the $x$ axes). In short $L(V)=2$ and since
$\mathcal{V}$ does not contain periodic orbits, the upper bound of two limit
cycles follows from Theorem 1.2.

Finally, we prove that one is the actual upper bound for the number of limit
cycles of (13). By Theorem 1.2 the stability of the limit cycles is given by
the sign of $-V(x,y)M_{1}(x)$ that coincides with the sign of $-H(x)>0.$ Hence
all of them are repelling hyperbolic limit cycles. By using the
Poincaré–Bendixson theorem it can be seen that the only situations compatible
with these results are that either (13) has no limit cycle, or that it has
exactly one, as we wanted to show.

The same tools allow to prove that

$\begin{cases}\dot{x}=y-bx^{3}(x-a)(2x-3a),\\\
\dot{y}=-x(x-a)(2x-a),\end{cases}$

has at most one limit cycle. In fact, notice that when $b=0$ it can be easily
integrated. It has two centers and a saddle point, and the separatrices of
this saddle point form two homoclinic trajectories which, together with the
critical point, have an eight shape. Numerically, the limit cycle seems to
bifurcate for $b\approx 0$ from this double loop and whenever it exists, it
surrounds the three equilibrium points of the system, a saddle, a strong focus
and a weak focus, both with different stabilities.

Similar examples to system (13), with more equilibrium points surrounded by a
limit cycle and for which Theorem 3.4 also works are not difficult to be
constructed. For instance, if we take $m=0,$ $F(x)=cx^{3}(1-x)(2-x)^{3}$ and
$G(x)=x^{2}(1-x)^{2}(2-x)^{2},$ with $|c|<2$, we get a system with five
critical points, two saddles and three foci, for which
$H(x)=8cx^{4}(1-x)^{4}(2-x)^{4}.$

Also, because it contains the classical van der Pol differential equation, we
particularize in detail a subfamily of the one given in item (ii). If we
consider $p=q=0,$ $c=1,$ and $A^{\prime}=C$ in $(ii)$ we get the following
result.

######  Corollary 3.5.

Consider the $\mathcal{C}^{1}$ differential system

$\begin{cases}\dot{x}=y-C(x)B(x),\\\ \dot{y}=-C(x)C^{\prime}(x),\end{cases}$

with $C(0)=0$ and $C^{\prime}(x)\neq 0$ for $x\neq 0.$ If $C(x)B^{\prime}(x)$
does not change sign and vanishes at isolated points, then this system has at
most one limit cycle and when it exists it is hyperbolic.

Notice that the van der Pol equation corresponds to $C(x)=x$ and
$B(x)=\lambda(x^{2}/3-1).$ Then $C(x)B^{\prime}(x)=2\lambda x^{2}/2$, which
does not change sign.

###### Proof of Corollary 3.5.

For these particular cases of differential systems contained in the family
(ii) we get that

$H(x)=2C^{3}(x)B^{\prime}(x).$

Hence, it does not change sign and vanishes at isolated points. Notice that
$C$ only vanishes at $x=0,$ because if $C(z)=0,$ by Rolle’s theorem
$C^{\prime}$ would vanish at a point between $0$ and $z.$ Hence,
$G^{\prime}(x)=2C(x)C^{\prime}(x)=0$ only at $x=0,$ and the corollary follows.
Observe also that in this case $V(x,y)=C^{2}(x)+y^{2}-yC(x)B(x),$ and the set
$\\{V(x,y)=0\\}$ has only one bounded connected component, the origin, and
then $L_{X}(V)=1.$ ∎

######  Remark 3.6.

For completeness we reproduce a second easy proof of the uniqueness and
hyperbolicity of the limit cycle of the van der Pol equation attributed to
Cherkas in [7, p. 105]. Write the equation as the system

$\begin{cases}\dot{x}=y,\\\ \dot{y}=-x-\lambda(x^{2}-1)y.\end{cases}$

By applying the Bendixson–Dulac theorem with $V=x^{2}+y^{2}-1$ and $s=2$ we
get that $M_{2}=2\lambda(x^{2}-1)^{2}.$ Clearly, the unit circle is not a
periodic orbit of the system, and $\\{V(x,y)=0\\}$ has two connected
components, one bounded and simply connected and a second one, say
$\mathcal{U},$ with $\ell(\mathcal{U})=1.$ Hence $L_{X}(V)=1$ and the result
follows.

When $b=0,$ the system introduced in item (iv) corresponds to the Wilson
Liénard equation ([28]) and when $|a|<2$ it has the algebraic limit cycle

$y^{2}-\frac{a}{6}x^{3}y+\frac{1}{144}(a^{2}x^{6}+144x^{2}-576)=0.$

Since this limit cycle is also hyperbolic we get that for $|b|$ small enough
the limit cycle persists and our theorem applies to get an upper bound of the
total number of limit cycles of the system when $b\neq 0$. We skip more
details because the study of this system is quite similar to the one that we
did for system (13).

We end this section studying in more detail the particular family of
differential systems of the form (11), introduced in [26],

$\begin{cases}\dot{x}=y-\lambda|y|^{m}(x^{3}-x),\\\ \dot{y}=-x,\end{cases}$
(14)

where $m\geq 2$ is an integer number and $\lambda\in\mathbb{R}.$ Notice that
the factor $x^{3}-x$ in (14) can be changed by $c^{2}x^{3}-x,$ with another
value of $\lambda,$ obtaining the same phase portrait. This is so, because by
doing the change of variables $(x,y)\to(cx,cy),$ with $c>0,$ the first
equation writes as $\dot{x}=y-c^{m}\lambda|y|^{m}(c^{2}x^{3}-x)$ and the
second one remains invariant. We do not take the factor as $x^{3}/3-x,$ which
corresponds to the van der the Pol equation when $m=0,$ simply to keep the
notation of [26]. We prove:

######  Theorem 3.7.

Consider the differential system (14) with $m\in\mathbb{N},$ and $m\geq 2.$
Then, the following holds:

  1. (i)

For $|\lambda|\neq 0$ small enough it has at least one limit cycle.

  2. (ii)

For $|\lambda|\geq\dfrac{3}{\sqrt{2}}\Big{(}\dfrac{3}{m}\Big{)}^{m/2}$ it has
no limit cycle.

###### Proof.

Notice that the case $\lambda=0$ corresponds to a linear center, and the phase
portrait when $\lambda<0$ can be easily obtained from the one with
$\lambda>0,$ simply by doing the change of variables and time
$(x,y,t)\to(x,-y,-t).$ Then, it suffices to make the proof for the case
$\lambda>0$.

$(i)$ Given any $\mathcal{C}^{1}$ perturbed Hamiltonian systems,

$\begin{cases}\dot{x}=\phantom{-}\dfrac{\partial H(x,y)}{\partial
y}+\varepsilon R(x,y),\\\\[10.0pt] \dot{y}=-\dfrac{\partial H(x,y)}{\partial
x}+\varepsilon S(x,y),\end{cases}$ (15)

where $\varepsilon$ is a small parameter, its associated
Melnikov–Poincaré–Pontryagin function is

$M(h)=\int_{\gamma(h)}S(x,y)\,{\rm d}x-R(x,y)\,{\rm d}y,$

where the curves $\gamma(h),$ for $h\in(h_{0},h_{1}),$ form a continuum of
ovals contained in $\\{H(x,y)=h\\}.$ It is known that each simple zero
$\bar{h}\in(h_{0},h_{1})$ of $M$ gives rise to a limit cycle of (15) that
tends, when $\epsilon\to 0,$ to $\gamma(\bar{h}),$ see for instance [8, 11].

Consider the differential system (14) with $\lambda=\varepsilon.$ By applying
the above result with $H(x,y)=x^{2}+y^{2}=h=r^{2},$ with $r\in(0,\infty),$ and
taking the parameterization of the level sets as $x=r\cos\theta,$
$y=r\sin\theta,$ we get that

$\displaystyle M(r^{2})$
$\displaystyle=\int_{x^{2}+y^{2}=r^{2}}|y|^{m}(x^{3}-x)\,{\rm d}y$
$\displaystyle=\int_{0}^{2\pi}\\!\\!\\!r^{m}|\sin\theta|^{m}(r^{4}\cos^{4}\theta-r^{2}\cos^{2}\theta)\,{\rm
d}\theta$
$\displaystyle=\frac{\sqrt{\pi}}{2}\frac{\Gamma\big{(}{(m+1)}/{2}\big{)}}{\Gamma\big{(}{(m+6)}/{2}\big{)}}r^{m+2}\big{(}3r^{2}-(m+4)\big{)},$

where $\Gamma$ is the Euler Gamma function. Hence, for each $m,$ this function
has a simple positive zero $r=\sqrt{{(m+4)}/{3}}$, that gives rise to the
desired limit cycle.

$(ii)$ We will apply Theorem 1.2 with $s=1/3$ and

$V(x,y)=\exp\left(\frac{\lambda^{2}y^{2m}}{9m}\right)\left(3+\lambda
xy|y|^{m-2}\right).$

Some calculations give that

$M_{1/3}=-\frac{1}{9}\exp\left(\frac{\lambda^{2}y^{2m}}{9m}\right){\lambda
x^{2}|y|^{m-2}}\left(2\lambda^{2}y^{2m}-27y^{2}+9(m-1)\right).$

We need that $M_{1/3}$ does not change sign. Hence, writing $y^{2}=z$ we want
that

$z^{m}-\frac{27}{2\lambda^{2}}z+\frac{9(m-1)}{2\lambda^{2}}\geq
0\quad\mbox{for}\quad z\geq 0.$ (16)

Let us prove, that given a real polynomial $P(z)=z^{m}+bz+c,$ with $m\geq 2,$
it holds that $P(z)\geq 0$ for all $z\geq 0$ if and only if

$b\geq-m\Big{(}\frac{c}{m-1}\Big{)}^{(m-1)/m}.$ (17)

Since $P(0)=c,$ an obvious first condition is that $c\geq 0.$ When $b\geq 0$
the result is trivial. When $b<0,$ since $p^{\prime}(z)=mz^{m-1}+b,$ the
function $P$ has a minimum at $z=z_{0}=(-b/m)^{1/(m-1)}.$ By imposing that
$P(z_{0})\geq 0,$ (17) follows after some straightforward computations.

Condition (17) applied to the polynomial (16) gives that

$-\frac{27}{2\lambda^{2}}\geq-m\left(\frac{9}{2\lambda^{2}}\right)^{(m-1)/m}.$

After some manipulations we get that this inequality is equivalent to the one
given in the statement.

Hence, we are under the hypotheses of Theorem 1.2. Moreover, since
$\\{V(x,y)=0\\}$ does not contain ovals, and all the connected components of
$\mathbb{R}^{2}\setminus\\{V(x,y)=0\\}$ are simply connected, we have that
$L_{X}(V)=0$ and the system has no limit cycle, as we wanted to prove. ∎

Similar computations that the ones of item (i) are done in the Appendix of
[26] for the case $m=2$.

######  Remark 3.8.

The result of item (ii) of Theorem 3.7 shows that for any $m\geq 2$ there
exits a value $\lambda=\lambda^{*}(m)$ such that for
$|\lambda|\geq\lambda^{*}(m)$ system (14) has no limit cycle. Moreover, it
gives an upper bound of this value.

For $m=2$ it is not a sharp bound, because in [26] the authors study
numerically the system and they find that $\lambda^{*}(2)\in(1.474,1.475),$
while our bound is $9\sqrt{2}/4\approx 3.182.$ Nevertheless, by using it, for
$m$ big enough, we prove that the limit cycles only exist for $\lambda$ in an
extremely thin interval of length $3\sqrt{2}\Big{(}\dfrac{3}{m}\Big{)}^{m/2}$
that decreases exponentially when $m$ grows.
It is the first time that the authors see a proof of the existence of this
type of exponentially small intervals for the presence of limit cycles.

###  3.3. About Massera’s theorem

Consider the classical Liénard equation

$\begin{cases}\dot{x}=y-F(x),\\\ \dot{y}=-x,\end{cases}$ (18)

with $F$ a class $\mathcal{C}^{2}$ function satisfying $F(0)=0.$ We prove, in
a very simple way, the following extension of the classical Massera’s theorem
([24, 25]), where the hyperbolicity of the limit cycle is also guaranteed.
Other authors had already proved this hyperbolicity, see for instance [19].

######  Theorem 3.9.

Consider the differential system (18). If the function $xF^{\prime\prime}(x)$
does not change sign and vanishes at isolated points, then it has at most one
limit cycle and when it exists it is hyperbolic.

###### Proof.

We will apply Theorem 1.2 with $V$ given by the function $K,$ defined in (5),
associated to the curvature of the system, and $s=1.$ By using the results of
Theorem 2.3 when $P=y-F(x)$ and $Q=-x,$ we obtain that

$V=K=x^{2}+y^{2}+F^{2}-2yF+x\big{(}y-F\big{)}F^{\prime}$

and $M_{1}=(y-F)^{2}xF^{\prime\prime}.$ Hence $M_{1}$ satisfies the hypothesis
of Theorem 1.2. To end the proof we have to show that $L_{X}(V)\leq 1.$ Since
the set $\mathcal{V}=\\{V(x,y)=0\\}$ does not contain orbits of the system, it
suffices to prove that $\mathcal{V}$ has at most one bounded connected
component, see Section 2.2. Clearly the points of $\mathcal{V}$ lie on the two
curves

$y=F(x)-\frac{1}{2}xF^{\prime}(x)\pm\frac{1}{2}\sqrt{x^{2}\big{(}(F^{\prime}(x))^{2}-4\big{)}}.$

Therefore the bounded connected components of $\mathcal{V}$ are given by $x=0$
and the bounded subsets of $\mathbb{R},$ where $(F^{\prime}(x))^{2}-4\geq 0.$
These components are either positively or negatively invariante by the flow of
the system because $M_{1}\big{|}_{\mathcal{V}}=\dot{V}$ does not change sign.
Hence they must surround some of the equilibrium points of the system. Since
the origin is the only equilibrium point, there is at most one of these
components. Hence, $L_{X}(V)\leq 1,$ as we wanted to prove. ∎

We want to emphasize the surprising simplicity of the proof of this classical
theorem with the methods employed in this work.

##  4\. Rigid systems

These systems write as

$\begin{cases}\dot{x}=-y+xF(x,y),\\\ \dot{y}=\phantom{-}x+yF(x,y),\end{cases}$
(19)

where $F$ is an arbitrary smooth function. This name is due to the fact that
in the usual polar coordinates $(r,\theta)$ it holds that $\dot{\theta}=1$
and, therefore, their flow rotates around the origin with constant angular
velocity, as a rigid rotation. Despite their simplicity and the fact that they
have the origin as the unique equilibrium point, the control of the number of
limit cycles of these systems is far to be completely known. They were
introduced by Conti in [2] and studied by several authors. We prove the
following result for them:

######  Theorem 4.1.

Let $X$ be the vector field associated to (19). If $F$ is of class
$\mathcal{C}^{2}$ and it holds that

$H:=F_{xx}F_{yy}-F^{2}_{xy}\geq 0,$ (20)

and $H$ vanishes on a null measure set, then (19) has at most $L_{X}(V)$ limit
cycles, where

$V=(x^{2}+y^{2})\left(xFF_{x}+yFF_{y}+xF_{y}-yF_{x}-1-F^{2}\right)$ (21)

and $L_{X}(V)$ is defined in Theorem 1.2.

###### Proof.

We apply again Theorem 1.2 with $V=K$ and $s=1,$ where $K$ is given in (5). We
can use the results of Theorem 2.3 with $P=-y+xF$ and $Q=x+yF.$ We get that
$V$ is as in (21) and

$M_{1}=D=(x^{2}+y^{2})\Big{(}\big{(}x^{2}F_{xx}+2xyF_{xy}+y^{2}F_{yy}\big{)}F^{2}\\\
+2\big{(}(x^{2}-y^{2})F_{xy}+xy(F_{yy}-F_{xx})\big{)}F\\\
+\big{(}x^{2}F_{yy}-2xyF_{xy}+y^{2}F_{xx}\big{)}\Big{)}.$

To control the sign of $M_{1}$ we first remove the factor $x^{2}+y^{2}.$
Notice that the discriminant of the remaining part, thinking it as a second
degree polynomial in $F$, $AF^{2}+BF+C,$ is
$B^{2}-4AC=-4(x^{2}+y^{2})^{2}H\leq 0.$ Moreover, looking to $A$ and $B$ as
quadratic homogenous polynomials of the form $ax^{2}+bxy+cy^{2},$ we get that
their corresponding discriminants coincide and are given by $b^{2}-4ac=-4H\leq
0.$ Therefore, the condition (20) implies that $M_{1}$ does not change sign
and vanishes only on a null measure set and hence our result follows. ∎

Notice that the upper bound for the number of limit cycles given in the above
theorem essentially depends on the shape of the set $\\{V(x,y)=0\\}.$ To get
the actual value of $L_{X}(V)$ for each case this set must be carefully
studied. We present now a concrete application when $F$ is a quadratic
polynomial.

######  Corollary 4.2.

Consider the rigid cubic system (19), with $F=a+bx+cy+dx^{2}+exy+hy^{2}.$ If
$4dh-e^{2}>0$ this system has at most one limit cycle, and when it exists it
is hyperbolic.

This result is not new. It was proved in [20] by using a totally different
approach: the authors transform the system into a periodic Abel differential
equation and then they apply know results about these equations. In that work
it is also proved that when $4dh-e^{2}<0$ there are systems with at least two
limit cycles. Our proof is different and self-contained. Another proof, based
on the study of the stability of the possible periodic orbits, is given in
[13].

###### Proof.

The function $H$ of Theorem 4.1 is $H(x,y)\equiv 4dh-e^{2}>0$ and hence the
system has at most $L_{X}(V)$ limit cycles. Here

$V(x,y)=(x^{2}+y^{2})\Big{(}-1-a^{2}+(c-ab)x-(ac+b)y+e(x^{2}-y^{2})+2(h-d)xy\\\
+(bx+cy)(dx^{2}+exy+hy^{2})+(dx^{2}+exy+hy^{2})^{2}\Big{)}.$

It is easy to verify that the set $\\{V(x,y)=0\\}$ does no contain orbits of
the system.

As the origin is the unique finite critical point of the system and $M_{1}$
does not change sign, the bounded connected components of the set
$\mathcal{V}=\\{V(x,y)=0\\}$ must surround the equilibrium point. In
principle, from the degree of $V(x,y)$ we can conclude that the maximum number
of them is three, being the origin one of these components. But it is easy to
show that there are at most two bounded connected components. This is so,
because if we take $y=0$ in the second factor of $V$ we obtain a polynomial in
$x$ of degree four, where the coefficient of $x^{4}$ and the independent term
are of opposite sign. Then it is not possible to have two positive roots and
two negative roots at the same time. Therefore, the number of connected
components in the set $\mathcal{V}$ is at most two, one of them being the
origin.

In the case where the second bounded connected component exists it is not
difficult to show that a limit cycle exterior to it cannot exist. The first
step is to determine the stability of infinity. Writing the system in polar
coordinates it is possible to show that, if $4dh-e^{2}>0$, the infinity is an
attractor for $d>0$ and it is repulsive for $d<0$. Moreover, it can be seen by
using that $M_{1}\big{|}_{\mathcal{V}}=\dot{V},$ that the flow associated to
the system traverses this second bounded component forward for $d>0$ and
inward for $d<0$. As between this bounded connected component of $\mathcal{V}$
and infinity only one limit cycle can exist, and if it exists it is
hyperbolic, taking into account the stability of infinity we conclude that no
limit cycle exists in this region. We conclude then that the system can have
at most one limit cycle. As the origin is an attractor for $a<0$ and it is
repulsive for $a>0,$ the limit cycle appears via a Hopf bifurcation at the
origin and must be located in the interior of the non trivial bounded
connected component of the set $\mathcal{V}.$ In conclusion, the limit cycle
exits if and only if $ad<0$ and it is unique. ∎

We end this work with a second corollary of Theorem 4.1 that also covers some
non-polynomial rigid differential systems.

######  Corollary 4.3.

Consider the differential system (19), with $F(x,y)=f(x)+g(y),$ where
$f(x)=\sum_{k=0}^{2n}f_{k}x^{k},$ with $f_{k}\geq 0,k\geq 2$ and $f_{2n}>0$.
We assume that $f^{\prime\prime}(x)\geq 0$, $g$ is of class $\mathcal{C}^{2},$
with $g^{\prime\prime}(y)\geq 0$ and it vanishes only at isolated points. Then
the system has at most two limit cycles. Moreover, if there exists $R>0$ such
that for all $(x,y)\in\mathbb{R}^{2}$ with $x^{2}+y^{2}\geq R^{2}$ it holds
that $F(x,y)\geq c>0,$ then the corresponding differential system has at most
one limit cycle, and when it exists, it is hyperbolic.

###### Proof.

For this case the function $H$ given in Theorem 4.1 is
$H(x,y)=f^{\prime\prime}(x)g^{\prime\prime}(y)\geq 0$ and it vanishes on a
null measure set. Hence the system has at most $L_{X}(V)$ limit cycles, where
$V$ is the function given in (21). To study $L_{X}(V),$ notice first that it
is not restrictive to assume that $g(0)=0.$ Then, $V(x,0)=x^{2}W(x),$ where

$W(x)=xf(x)f^{\prime}(x)+g^{\prime}(0)x-1-f^{2}(x).$

Since when $f(x)=f_{k}x^{k}$ it holds that
$xf^{\prime}(x)-f(x)=(k-1)f_{k}x^{k},$ we get easily that
$W(x)=\sum_{j=0}^{4n}w_{k}x^{k},$ where all $w_{k}\geq 0$ for $k\geq 2$ and
$w_{0}=-1-f^{2}(0)<0.$ Hence, by the Descarte’s rule of signs the number of
positive roots of $W$ is 1. As a consequence, the set
$\mathcal{V}=\\{V(x,y)=0\\}$ has at most one bounded connected component
surrounding the origin, different from the origin itself. Recall that $V$ has
the factor $x^{2}+y^{2}.$ Notice that the set $\mathcal{V},$ which is not
invariant by the flow of $X,$ can not contain other bounded connected
components. This is so, because $M_{1}$ does not change sign, and
$M_{1}\big{|}_{\mathcal{V}}=\dot{V}.$ Therefore, these connected components
must surround some equilibrium point of $X,$ but the origin is the only one.
As a consequence of the above reasoning $L_{X}(V)\leq 2,$ see Section 2.2.

Let us prove now that under the hypothesis on the growth of $F$ the maximum
number of limit cycles is 1. Notice that if $r=\sqrt{x^{2}+y^{2}},$ it holds
that $\dot{r}=rF(r\cos\theta,r\sin\theta)\geq cr.$ Hence, the infinity is an
attractor, and we can use similar arguments that in the proof of Corollary 4.2
to show that when $L_{X}(V)=2$ the differential system has no limit cycle in
the unbounded components of $\mathbb{R}^{2}\setminus\mathcal{V}.$ Therefore,
the existence of at most one limit cycle follows. ∎

### Acknowledgments

This work has received funding from the Ministerio de Ciencia e Innovación
(PID2019-104658GB-I00 grant) and the Agència de Gestió d’Ajuts Universitaris i
de Recerca (2017 SGR 1617 grant).

## References

  * [1] Bendixson, I. Sur les courbes définies par des équations différentielles.  Acta Math. 24, (1901) 1–88.
  * [2] Conti, R. Uniformly isochronous centers of polynomial systems in ${\bf R}^{2}$.  In Differential equations, dynamical systems, and control science, vol. 152 of Lecture Notes in Pure and Appl. Math. Dekker, New York, 1994, pp. 21–31.
  * [3] Chen, H., and Chen, H. Global dynamics of a Wilson polynomial Liénard equation.  Proc. Amer. Math. Soc. 148, (2020) 4769–4780.
  * [4] Cherkas, L.A. Dulac function for polynomial autonomous systems on a plane.  Differ. Equ. 33, (1997) 692–701.
  * [5] Cherkas, L.A., and Grin’, A.A. On the Dulac function for the Kukles system.  Differ. Equ. 46, (2010) 818–826.
  * [6] Chicone, C. Bifurcations of nonlinear oscillations and frequency entrainment near resonance.  SIAM J. Math. Anal. 23, (1992) 1577–1608.
  * [7] Chicone, C. Ordinary differential equations with applications. In: Texts in Applied Mathematics, vol. 34, 2nd edn. Springer, New York (2006).
  * [8] Christopher, C., and Li, C. Limit cycles of differential equations.  Advanced Courses in Mathematics. CRM Barcelona. Birkhäuser Verlag, Basel, 2007.
  * [9] Coppel, W. A. A survey of quadratic systems.  J. Differential Equations 2, (1966) 293–304.
  * [10] Diliberto, S. P. On systems of ordinary differential equations.  Contributions to the Theory of Nonlinear Oscillations, pp. 1–38.  Annals of Mathematics Studies, no. 20. Princeton University Press, Princeton, N.J., 1950.
  * [11] Dumortier, F., Llibre, J., and Artés, J. C. Qualitative theory of planar differential systems.  Universitext. Springer-Verlag, Berlin, 2006.
  * [12] Dulac, H. Récherches des cycles limites.  C.R. Acad. Sci. Paris Sér. I Math. 204, (1937) 1703–1706.
  * [13] Freire, E., Gasull, A., and Guillamon, A. Limit cycles and Lie symmetries.  Bull. Sci. Math. 131, (2007) 501–517.
  * [14] Garcia, R.A., Gasull, A., and Guillamon, A. Geometrical conditions for the stability of orbits in planar systems.  Math. Proc. Cambridge Philos. Soc. 120, (1996) 499–519.
  * [15] Gasull, A., and Giacomini, H. A new criterion for controlling the number of limit cycles of some generalized Liénard equations.  J. Differ. Equations 185, (2002) 54–73.
  * [16] Gasull, A., and Giacomini, H. Upper bounds for the number of limit cycles through linear differential equations.  Pacific J. Math. 226, (2006) 277–296.
  * [17] Gasull, A., and Giacomini, H. Upper bounds for the number of limit cycles of some planar polynomial differential systems.  Discrete Contin. Dyn. Syst. 27, (2010) 217-229.
  * [18] Gasull, A., and Giacomini, H. Some applications of the extended Bendixson–Dulac theorem.  In Progress and challenges in dynamical systems, vol. 54 of Springer Proc. Math. Stat. Springer, Heidelberg, 2013, pp. 233–252.
  * [19] Gasull, A., Giacomini, H., and Llibre J. New criteria for the existence and non-existence of limit cycles in Liénard differential systems.  Dyn. Syst. 24, (2009) 171–185.
  * [20] Gasull, A., Prohens, R., and Torregrosa, J. Limit cycles for rigid cubic systems.  J. Math. Anal. Appl. 303, (2005) 391–404.
  * [21] Gasull, A., and Sabatini, M. Fixed and moving limit cycles for Liénard equations.  Ann. Mat. Pur. Appl. 198, (2019) 1985–2006.
  * [22] Ilyashenko, Yu. Centennial history of Hilbert’s 16th problem.  Bull. Am. Math. Soc. (N.S.) 39, (2002) 301–354.
  * [23] Lloyd, N. G. A note on the number of limit cycles in certain two-dimensional systems.  J. London Math. Soc. (2) 20, (1979), 277–286.
  * [24] Massera, J. L. Sur un théorème de G. Sansone sur l’équation di Liénard.  Boll. Un. Mat. Ital. 9, (1954) 367–369.
  * [25] Sansone, G. Soluzioni periodiche dell’equazione di Liénard. Calcolo del periodo.  Rend. Sem. Mat. Univ. Politec. Torino 10, (1951), 155–171.
  * [26] Villari, G, and Zanolin, F. On the qualitative behavior of a class of generalized Liénard planar systems. Preprint 2020\.
  * [27] Wilson, G. Hilbert’s sixteenth problem.  Topology 17, (1978) 53–73.
  * [28] Wilson, J. C. Algebraic periodic solutions of $\ddot{x}+f(x)\dot{x}+g(x)=0.$ Contrib. Differ. Equ. 3, (1964) 1–20.
  * [29] Yamato, K. An effective method of counting the number of limit cycles.  Nagoya Math. J. 76, (1979) 35–114.
  * [30] Ye, Yan Qian; et alt. Second edition. Translations of Mathematical Monographs, 66. American Mathematical Society, Providence, RI, 1986.