Guidelines for benchmarking of optimization approaches for fitting
mathematical models

Clemens Kreutz1,2*

1 Institute of Medical Biometry and Statistics, Faculty of Medicine and
Medical Center, University of Freiburg, Stefan-Meier-Str. 26, 79104 Freiburg,
Germany
2 CIBSS – Centre for Integrative Biological Signalling Studies, University of
Freiburg, 79104 Freiburg, Germany

*<EMAIL_ADDRESS>
## Abstract

Insufficient performance of optimization approaches for fitting of
mathematical models is still a major bottleneck in systems biology. In this
manuscript, the reasons and methodological challenges are summarized as well
as their impact in benchmark studies. Important aspects for increasing
evidence of outcomes of benchmark analyses are discussed. Based on general
guidelines for benchmarking in computational biology, a collection of tailored
guidelines is presented for performing informative and unbiased benchmarking
of optimization-based fitting approaches. Comprehensive benchmark studies
based on these recommendations are urgently required for establishing of a
robust and reliable methodology for the systems biology community.

## Introduction

A broad range of mathematical models are applied in systems biology. Depending
on the questions of interest and on the amount of available data, the type of
models and the level of detail vary. Most frequently, _ordinary differential
equation models (ODEs)_ are applied because they enable a non-discretized
description of the dynamics of a system and allow for quantitative evaluation
of experimental data including statistical interpretations in terms of
confidence and significance. In the BioModels Database [22], currently 83% of
all models which are uniquely assigned to a modelling approach are ODE models.
In this manuscript, we focus on optimization-based fitting of these models
although many aspects are general, and thus also relevant for other modelling
types and approaches.

Abundances of compounds, the strengths and velocities of biochemical
interactions are typically context-dependent, i.e. vary between species,
tissues and cell types. Thus, they are represented as unknown parameters in
mathematical models. Application-specific calibration of these models is
therefore required which corresponds to estimation of these unknown parameters
based on experimental data.

In most cases, parameter estimation is performed by optimization of a suitable
_objective function_ like minimization of the sum of squared residuals for
_least-squares_ estimation, or maximization of the likelihood for _maximum
likelihood estimation_ [2]. Although parameter estimation is a central task of
modelling, lack of reliable computational approaches for fitting is still a
bottleneck in systems biology. The absence of high-performing software
implementations seems a major reason, why ODE-based modelling is not yet a
routinely applied computational approach for analyzing experimental data.

The importance of proper designs for benchmark studies in computational
biology has been discussed in several publications [5, 19, 20]. General
guidelines have been provided recently for computational analysis of omics
data [23], for multiple alignment of protein sequences [1], for supervised
classification methods [6] as well as for periodic scientific benchmarking [7]
and for general studies in computational biology [4, 27, 36]. In this
manuscript, these aspects are discussed in the context of benchmarking
optimization approaches for fitting mathematical models in systems biology.

## Why is reliable fitting challenging?

A major characteristic of the mathematical models applied in systems biology
is the intention to mirror the biological process of interest because this
facilitates enhanced possibilities of interpretations and understanding. For
this purpose, molecular compounds like proteins and spatial compartments like
cells are defined as model components and regulatory interactions between the
considered compounds are translated as interactions into the model. In
contrast to _phenomenological_ models which describe empirical relationships
in an abstract and simplified manner, the complexity of these so-called
_mechanistic_ models is dictated by the complexity of the investigated
biological process.

In systems biology, the amount of available experimental data for parameter
estimation is always limited. In such settings, multiple parameter
combinations can yield to the same model response for experimentally
investigated conditions. A measured steady-state concentration, as an example,
might only provide information about the ratio $k_{prod}/k_{deg}$ of
production and degradation rates. Because, different combinations of the
individual parameters $k_{prod}$ and $k_{deg}$ yield the same steady state,
all combinations with the same ratio fit the data equally well. This
characteristic has been termed _non-identifiability_. Because of non-
identifiability, ill-conditioning and entirely flat manifolds are common
features of systems biological models.

The following typical attributes of mechanistic models in systems biology
raise methodological challenges for model fitting:

  1. 1.

The models are large in terms of number of parameters. Thus, optimization has
to be performed in high-dimensional spaces.

  2. 2.

The objective functions used for fitting depend (strongly) non-linear on the
model parameters. Therefore, numerical optimization is intricate and local
optima can exist.

  3. 3.

Evaluation of the objective function comprises numerical integration of the
ODEs which is computationally demanding and only feasible with limited
numerical accuracy.

  4. 4.

Integration of the ODEs cannot be performed analytically. Therefore, exact
mathematical calculations are infeasible.

  5. 5.

Derivatives of the objective function have limited numerical accuracy and
cannot be calculated naively.

  6. 6.

Parameter values vary over several orders of magnitudes and usually only a
limited amount of prior knowledge is available. Therefore, it is difficult to
specify priors, initial guesses, and/or bounds.

  7. 7.

Optimization has to cope with constraints like upper and lower bounds and with
non-identifiability, i.e. with ill-conditioning and entirely flat manifolds.

In addition, discontinuities of external inputs (so-called _events_) [15]
might occur which has to be handled properly. Sometimes, steady state initial
values have to be calculated numerically which is another source of numerical
inaccuracy [12].

## Existing approaches and benchmark studies

Figure 1: Tasks to be accomplished for fitting ODE models. Fitting of ODE
models requires several tasks. Most approaches combine both, global and local
search strategies. A prominent global search strategy is random drawing of
multiple initial guesses and performing local optimization for each starting
point. In each optimization step, the ODEs have to be solved for evaluation of
the objective function $\chi^{2}(\theta)$. Moreover, an incremental
improvement strategy is required for suggesting a new parameter vector for the
next iteration step. This task is usually performed based on derivatives or
approximations thereof.

For fitting of ODE models several computational tasks have to be combined as
indicated in Figure 1. For these tasks, a multitude of methods have been
published (e.g. as summarized in [2, 14, 33]). For a given parameter vector,
the ODEs have to be solved by numerical ODE integration methods in order to
evaluate the objective function. Moreover, iterative optimization requires a
strategy for suggesting the next trial parameters in each optimization step.
This requires calculation or approximation of derivatives or alternative
approaches for incremental improvement. Such iterative local optimization
approaches are typically deterministic and are usually combined with a
stochastic global search strategy in order to enable convergence to the global
optimum.

Currently, the existing benchmark studies provide only a heterogeneous,
fragmentary and partly inconsistent picture about the performance and
applicability of fitting techniques and how these different numerical tasks
have to be efficiently combined.

An early publication in systems biology found superior performance of so-
called _multiple shooting_ [26] for local optimization. This method, however,
is difficult or even impossible to apply for partly observed systems with
complex observation functions. It requires a custom implementation of ODE
integration method and combination with global search strategies was not
considered. Moreover, there are no publicly available implementations.
Consequently, this methodology disappeared from the systems biology field
during recent years.

Multi-start deterministic optimization showed superior performance in the
_Dialogue for Reverse Engineering Assessments and Methods (DREAM)_ benchmark
challenges about parameter estimation [32] and network reconstruction [24].
Here, a _trust-region_ and gradient based deterministic nonlinear least
squares optimization approach [8] has been utilized as local optimization
strategy and global search was performed by utilizing multiple runs with
random initial guesses. This approach is implemented in the _Data2Dynamics_
modelling framework [30] and has been shown as superior to other approaches in
several studies [29, 37, 31].

In [29] it has been shown that deterministic gradient-based optimization is
superior to a batch of hybrid and pure stochastic algorithms. In contrast,
[10, 11, 16, 25] found superior performance of stochastic optimization
methods. In [34] it was confirmed that multi-start gradient-based local
optimization is often a successful strategy but, on average, a better
performance could be obtained with a hybrid metaheuristic combining
deterministic gradient-based optimization with a global scatter search
metaheuristic [10].

In [29, 13], it was shown that derivative calculation based on finite
difference is inappropriate. This outcome is partly questioned in [9], at
least if optimization is performed on a normalized data scale. [13] showed
that so-called _adjoint-sensitivities_ are computationally most efficient for
derivative calculations for large models.

It has been repeatedly claimed that parameters are preferably optimized on the
log-scale [20, 34, 17]. Nevertheless, optimization of model parameters is
still frequently performed at the linear scale in applications and even in
benchmark studies like [9].

In addition to neutral publications which focus primarily on benchmarking,
there are additional papers which focus on introducing a new approach and
compare the performance of alternative approaches in less comprehensive
manner. Despite this rather large amount of studies, there is currently no
consensus and there are no clear rules in the systems biology community about
proper selection of approaches for fitting. This shows that up to now,
benchmark studies do not provide convincing evidence which reveals the
requirement for improved benchmark analyses.

## Pitfalls of benchmark studies

### P1: Unrealistic setup

Realistic benchmark studies have to utilize restrict to the identical setting
and to the same amount of information which is available in an application
setting. Simulated data deviate from this ideal setting because real data in
molecular biology typically contain non-trivial correlations, artifacts or
systematic errors which is usually not considered for simulated data sets.
Moreover, for rejecting incomplete model structures, fitting of experimental
data must also work for wrong models. Since simulated data is typically
generated with the same model structure used for fitting, there is usually no
mismatch between model and data, and thus optimization is only evaluated for
settings where a correct model structure is available. In order to not rely on
such critical assumptions, it is strongly preferable to assess fitting
approaches based on real experimental data.

In real applications, there is no information available for tuning of
configuration parameters like integration tolerances or thresholds defining
termination criteria for iterative optimization. Therefore, default
configuration parameters have to be utilized or they have to be defined by
consistency checks. In order to prevent unrealistic performance assessment,
one has to strictly avoid tuning of configuration parameters based on
performance criteria which are not available in practice. Instead, it should
be prespecified before evaluation of the methods, how configuration parameters
will be chosen. Moreover, the process of tuning these parameter has to be
counted as additional runtime.

### P2: Ignoring cofactors

Abbrev. | Cofactor | Typical possible choices
---|---|---
C1 | Application problem | Model equations and the data set(s)
C2 | Primary performance criteria | Convergence per computation time, iteration steps
C3 | Secondary performance criteria | Documentation, user-friendlyness, code quality
C4 | Parameter scale | Linear vs. log-scale
C5 | Global search strategy | Multiple initial guesses, scatter search algorithms, …
C6 | Initial guess | Fixed vs. random
C7 | Initial guess distribution | Normal vs. uniform vs. latin-hypercube
C8 | Parameter constraints | Upper and lower bounds
C9 | Prior knowledge | None vs. (weakly) informative priors
C10 | ODE integration implementation | SUNDIALS, Matlab, R, …
C11 | ODE integration algorithm | Stiff vs. non-stiff approaches, Adams-Moulton vs. BDF, …
C12 | Integration accuracy | ODE integrator tolerances
C13 | Derivative calculation | Finite differences, sensitivity equations, adjoint sensitivities
C14 | Stopping rule | Optimization termination criteria
C15 | Handling of non-converging integration | Termination of optimization vs. infinite loss
C16 | Algorithm-specific configurations | Cross-over rate, annealing temperature, number of particles, …
Table 1: Cofactors. The performance of an optimization approaches depends on
many decisions and configurations C1-C16. For comparison of several
approaches, these attributes appear as cofactors. Performance benefits for
individual choices do not necessarily indicate a general advantage because
benefits might originate from the chosen configurations.  Figure 2: Impact of
configuration parameters. Multiple configurations can have an impact on the
performance and performance benefits. In this illustration example, two
optimization approaches have different sensitivity with respect to the choice
of tolerances controlling the numerical error of ODE integration. Moreover,
both approaches have different optimal choices for this tuning parameter. For
most integration tolerances, approach B is superior. However, the overall best
performance has approach A for optimally chosen tolerances. This illustration
example highlights the importance of evaluation of configuration parameters
for drawing valid conclusions.

A lot of decisions have to be made for benchmarking of optimization
approaches. First, application problems have to be selected. Then, a strategy
for combining global and local search has to be specified and an ODE
integration method has to be selected. Moreover, parameter bounds and
parameter scales (linear vs. logarithmic) have to be defined, ODE integration
algorithms and tolerances have to be specified as well as stopping criteria
for iterative optimization. Each of these decisions or _cofactors_ can have a
remarkable impact on the outcome of a benchmark study.

The severe problem is that these factors affect the performance of individual
optimization approaches. In the example illustrated in Figure 2, approach A
requires properly tuned tolerances controlling the accuracy of ODE
integration. If properly chosen, the approach is superior. In contrast,
approach B is less sensitive to these tolerances and outperforms A for most
choices although the approach can never reach the maximal performance which is
possible for approach A.

Thinking in a multi-factorial terminology is not yet common for benchmark
studies, but well-established e.g. for clinical studies where patient-specific
covariables like gender, age, smoking, etc. occur as cofactors. Table 1
summarizes 16 typical cofactors for optimization benchmark studies and
illustrates that many decisions have to be made to implement an optimization-
based parameter estimation approach.

The impact of cofactors can be minimized by proper study designs. Moreover,
one can estimate the impact of cofactors on optimization performance by
multivariate analysis of benchmarking outcomes. For the illustration example
in Figure 2, one could determine the impact of the cofactor by evaluation of
the whole range in the benchmark study which is also feasible in real
application settings. Treatment of cofactors in optimization benchmark studies
has been discussed in [20].

### P3: Performing only case studies

One of the most important cofactor is the chosen set of application problems,
i.e. the models and the data sets which are used for benchmarking. A chosen
model determines the intricacy of the optimization task in terms of dimension
(number of parameters), non-linearity, ill-conditioning, local optima, amount
of information provided by the data etc. and thereby affect performance.

If only a single model is evaluated, the study has only low evidence because
the behaviour of the applied approaches might be completely different for
another application problem. In terms of evidence, such a study is only a
_case report_. Therefore, benchmark studies should be performed based a
comprehensive set of test problems in order to draw generally valid
conclusions.

### P4: Non-convergence and local optima are hardly distinguishable

Figure 3: Ambiguous interpretation of optimization outcomes. For non-trivial
optimization problems, the results of independent optimization runs are
typically not the same. The upper left panel indicates an outcome for three
optimization runs, e.g. generated with different starting points. If the
objective function values after optimization are different (“scenario A”),
such an outcome could be explained by local optima (“explanation A1”) or by
convergence problems of the optimization algorithm. If the same values for
objective function are obtained, there might be several local optima with the
same objective function (“explanation B1”), there might be a convergence
problem (“B2”) or non-identifiabilities might exist (“B3”), i.e. the estimated
parameters are not uniquely specified by the data and then entire manifolds
yield the same objective function.

Sub-optimal convergence behavior is difficult to be discriminated from local
optima since in both cases optimization terminates at distinct points in the
parameter space, usually with different objective function values. In high-
dimensional spaces, it is difficult to evaluate whether a point in the
parameter space is a local optimum, especially if the objective function and
its derivatives can only be evaluated with limited numerical accuracy.
Existing approaches which could be applied like the profile likelihood [28],
reconstruction of flat manifolds [18], identifiability analysis [21], or
methods based on a local approximation of the Hessian do not enable reliable
classification in case of convergence issues and/or numerical inaccuracies.
Thus, convergence issues can easily be confused with local optima. New methods
for proving local optimality and for detection of convergence issues would be
very valuable for the assessment and for future improvements of optimization
approaches.

The upper left panel in Figure 3 shows a hypothetical outcome from three
optimization runs. “Scenario A” indicates that both, convergence problems and
local optima might generate the same observations. Even if identical values
for the objective function are obtained (“scenario B”), the interpretation is
ambiguous because there are no reliable approaches to distinguish convergence
issues from local optima. Moreover, a third possible explanation in this
scenario are non-identifiabilities, i.e. multi-dimensional optimal sub-spaces.
If non-identifiabilities exist, even high-performing optimization approaches
converge to distinct points in the parameter space. For benchmark analyses,
this means that measuring similarity of the outcomes of multiple optimization
runs by distances in the parameter space is not reasonable.

## Guidelines for benchmark studies

Ten general guidelines for benchmark studies in whole computational biology
field were recently presented in [36]. Here, we concretize these suggestions
for benchmarking of parameter optimization approaches and discuss their
relevance and feasibility. Our tailored guidelines G1-G7 and G9 correspond to
the recommendations presented in [36], guideline G8 comprises two distinct
guidelines (no. 8 and 10 in [36]). In order to provide readers a summary about
the design of a benchmark study and the resulting evidence in an easily
accessible and clearly structured manner, we recommend that the guidelines
presented in the following are discussed in future publications in a point-by-
point manner.

### G1: Clear definition of aim and scope

Benchmark studies which are performed to illustrate benefits of a newly
presented approach are easily biased because they are performed to confirm
merits. Often, application problems are utilized where existing methods have
limited performance. Consequently, it should be precisely stated, whether
benchmark analyses are performed for introducing a new approach, or whether
the goal is performing a neutral and comprehensive study based on previously
published computational approaches and test cases.

Like in other fields of computational biology, the performance of optimization
methods is context-specific, i.e. depends on the chosen benchmark problem.
Thus, it is essential to define the scope of a study and select representative
test cases according to this definition [20]. For mathematical modelling in
the systems biology field, one could define the scope by the biological
background. Models of signalling pathways, as an example, usually describe the
dynamics of activation after stimulation and show transient dynamic responses.
In contrast, metabolic models are usually steady state descriptions. Gene
regulatory networks, on the other hand, typically have distinct rate laws
because activatory and inhibitory effects are described by products of _Hill
equations_ , instead of _mass action kinetics_. Other possibilities for
defining scopes could be based on model attributes like amount of data, number
of parameters, existence of events or steady state constraints, or based on
dynamic characteristics like occurrence of oscillations.

### G2: Inclusion of (all) relevant methods

Published optimization methods are usually only available in distinct software
or programming environments and require combination of several numerical tasks
which do not perfectly coincide in different programming languages or software
package. As an example, the trust-region based nonlinear least-square
optimization approaches implemented in Matlab and R are not identically
programmed and produce different outcomes with differing performances [38].

Moreover, there is not yet an established standard for defining an
optimization problem comprehensively including all model equations, data,
measurement errors, priors, constraints, algorithms, and configuration
parameters. Thus, it is very elaborate and partly infeasible to include
approaches which are previously applied in other benchmark studies. This
limitation demands for standardized data and documentation formats like
_COMBINE_ archives [3].

Nevertheless, it should be aspired to include as many approaches as possible.
For reasonable interpretation of observed performances, it absolutely is
essential to implement at least one state-of-the-art approach.

### G3: Selection of representative test cases

Because the performance of optimization approaches strongly differ between
application models, a rather large number of models is required to obtain a
representative and comprehensive picture. However, the number of available
benchmark problems is up to now strongly limited. Six benchmark models have
been published in [35], however four of them only contain simulated data.
Recently, a set of 20 benchmark problems with experimental data sets has been
published [17], but it still remains difficult to perform larger studies.
Because of the small number of available test cases, it is currently only
hardly feasible to define a narrow scope while still having enough test
problems.

Simulating data for assessing modelling approaches requires much more
specifications than in most other fields of computational biology because
realistic combinations of sampling times, observables, observation functions,
error models, and experimental conditions have to be defined. Moreover, as
argued above as pitfall P1, simulation data only has limited value for
benchmarking of optimization approaches. Therefore, benchmark studies based on
experimental data are currently strongly recommended.

### G4: Appropriate configuration parameters and software versions

Assessment of fitting approaches requires many decisions at the level of
selecting approaches for the different numerical tasks as well as on the level
of configuration parameters. As discussed above (pitfall P2), configuration
parameters can have a major impact on performance. In order to guarantee fair
comparisons, the strategy for their definition has to be pre-specified and
should be done in an unbiased manner, i.e. it has to be carefully evaluated
that individual optimization approaches are not privileged.

According to our knowledge, existence and impact of differences between
subsequent versions of the same software package are not yet investigated.
Nevertheless, it is essential to comprehensively describe the applied approach
including software versions in order to guarantee full reproducibility.
Software tools for enhancing reproducibility of computational analyses were
summarized in [36].

### G5: Evaluation in terms of key quantitative metrics

Optimization is in almost all circumstances assessed by means of convergence,
i.e. in terms of probabilities or frequencies of finding local or global
optima. Although all local minima with statistically valid objective function
are of interest, commonly the primary interest is the global optimum. For
proving applicability and testing performance, it is usually also sufficient
to consider convergence to any kind of optimum because whether the global or a
local optimum is identified is often only a question of the chosen initial
guess and strongly depends on the size of the search space. For deterministic
optimization, as an example, each optimum has a region of attraction. Hence,
the frequency of finding the global optimum relative to the frequency of
converging to a local optimum is mainly a question of the size of the search
space and location of the optimization starting points.

Computational efforts have to be distributed among global and local search
strategies, e.g. computational efforts can be spend either for an increasing
number or for increasing lengths of the individual optimization runs. In order
to balance this trade-off, it is a very reasonable strategy to assess
convergence to a local/global optima by calculating the expected runtime for a
single converged run as recently suggested in [34]. It should be kept in mind,
however, that runtimes are also dictated by the computer system and especially
on parallelization. Thus, one has to ensure that rating of runtimes are fair.
Moreover, it should be investigated whether outcomes depend on the amount of
parallelization, i.e. on the number of processors, or on switching on/off
parallelized code execution.

An additional issue is the definition of convergence which is typically done
by thresholds for the objective function relative to the overall best known
solution. The choice of the threshold is a cofactor and its impact has to be
investigated. Moreover, the order of magnitude has to be chosen properly, i.e.
to guarantee that only fits which are in statistical agreement with data and
measurement errors are counted as successfully converged.

In most other benchmarking fields, one advantage of simulated data is
knowledge about the underlying truth because then appraisal is feasible in
terms of true/false positives or in terms of bias and variance. This is
different for benchmarking of optimization approaches. On the one hand, the
global solution might even for simulated data unknown because each noise
realization has a different optimum. On the other hand, outcomes from several
optimization approaches can be assessed by evaluating the objective function
which is feasible without any restriction for experimental data.

### G6: Evaluation of secondary measures

The quality of documentation, user-friendlyness, and code quality are
important for diminishing the risk of incorrect application of optimization
approaches. These attributes are therefore valuable as secondary measures.

Other secondary measures mentioned in the literature, are only subordinately
relevant since bad convergence behavior cannot be compensated by other
aspects. Moreover, traditional trade-offs like between precision and recall or
between bias and variance do not apply for assessing convergence of
optimization algorithms.

Feasibility of adaptation of optimization algorithms and fitting approaches by
users is also not a major secondary aim because custom heuristics are not
recommended from our perspective. Within a model calibration approach, there
are many aspects and details which comprehensively have to match together.
Therefore, algorithm development requires expert knowledge and traditional
strategies for tuning optimization approaches should be exploited until
sufficient convergence behavior is obtained instead of implementing a weakly
tested custom solution.

### G7: Interpretation and recommendation

Fitting of ODE models requires several numerical tasks. Optimization only
works reliably, if all these components match together and perform
sufficiently well. Proper interpretation thus requires evaluation of the
impact of all configuration options (see pitfall P2) and the multi-factorial
setup has to be accounted [20].

Benchmark studies should provide clear recommendation rules about the
selection of optimization approaches for specific scopes of application, e.g.
by deriving decision trees like “use approach A in case X, use B otherwise”.

An advantage of benchmark studies in this field is that multiple optimization
approaches can be applied subsequently in order to optimize the objective
function. One could therefore apply multiple well-performing approaches
consecutively if explicit rules for selecting single approaches cannot be
derived.

### G8: Publication and reproducible reporting of results

Since for optimization small technical details might deteriorate performances
and thereby strongly impact the outcomes of benchmark studies, it is
critically essential to publish the source code for the performed analyses as
well as details and versions of programming environments, operating systems,
and software packages.

In order to obtain consensus within the community, it seems very promising to
implement online resources for fitting. According to [36], we also recommend
that pre-print versions of peer-reviewed articles should be published in order
to speed up the distribution of results, broaden accessibility, and permitting
additional feedback.

### G9: Enabling future extensions

It is highly valuable to design and perform benchmark analyses which allow
extensions by additional optimization approaches as well as new test problems
in order to permit continuous update and refinements and to avoid that studies
become outdated. This guideline presented in [36] also fully applies to
benchmarking of optimization approaches.

## Conclusions

For optimization-based parameter fitting of ODE models, there are several
methodological challenges. According to existing benchmark studies, it seems
that ODE models from the systems biology field demand such a rich variety of
methodological requirements that each optimization approach can fail.
Unfortunately, existing studies provide only a fragmentary and inconsistent
picture about the applicability of existing approaches and there is no
consensus about proper selection of optimization approaches in the systems
biology field. Thus, reliable fitting of mathematical models remains a
limiting bottleneck in systems biology.

We presented pitfalls for the design, analysis and interpretation of benchmark
studies. Moreover, general guidelines from the literature were tailored to the
optimization-based parameter estimation setting. The presented pitfalls and
guidelines also indicate further methodological requirements. Standards for
exchanging model analyses have to be promoted in order to permit comparisons
over all software environments. A reliable methodology for distinguishing
numerical causing convergence problems from local optima has to be
established. Moreover, there is urgent need for more benchmark problems which
could be partly resolved by developing methodology for simulating data in a
realistic manner.

## Abbreviations

G: Guideline; ODEs: Ordinary differential equations; P: Pitfall

## Acknowledgements

This work was supported by the German Ministry of Education and Research by
grant EA:Sys _FKZ031L0080_ and by the German Research Foundation (DFG) under
Germany’s Excellence Strategy _CIBSS-EXC-2189-2100249960-390939984_.

## References

  * 1. M. R. Aniba, O. Poch, and J. D. Thompson.  Issues in bioinformatics benchmarking: the case study of multiple sequence alignment.  Nucleic Acids Research, 38(21):7353–7363, 2010.
  * 2. M. Ashyraliyev, Y. Fomekong-Nanfack, J. A. Kaandorp, and J. G. Blom.  Systems biology: parameter estimation for biochemical models.  The FEBS journal, 276(4):886–902, 2009.
  * 3. F. T. Bergmann, R. Adams, S. Moodie, J. Cooper, M. Glont, M. Golebiewski, M. Hucka, C. Laibe, A. K. Miller, D. P. Nickerson, et al.  Combine archive and omex format: one file to share all information to reproduce a modeling project.  BMC Bioinformatics, 15(1):369, 2014.
  * 4. A.-L. Boulesteix.  Ten simple rules for reducing overoptimistic reporting in methodological computational research, 2015.
  * 5. A.-L. Boulesteix, H. Binder, M. Abrahamowicz, W. Sauerbrei, and S. P. of the STRATOS Initiative.  On the necessity and design of studies comparing statistical methods.  Biometrical Journal, 60(1):216–218, 2018.
  * 6. A.-L. Boulesteix, S. Lauer, and M. J. Eugster.  A plea for neutral comparison studies in computational sciences.  PloS one, 8(4):e61562, 2013.
  * 7. S. Capella-Gutierrez, D. De La Iglesia, J. Haas, A. Lourenco, J. M. F. Gonzalez, D. Repchevsky, C. Dessimoz, T. Schwede, C. Notredame, J. L. Gelpi, et al.  Lessons learned: Recommendations for establishing critical periodic scientific benchmarking.  BioRxiv, page 181677, 2017.
  * 8. T. F. Coleman and Y. Li.  An interior trust region approach for nonlinear minimization subject to bounds.  SIAM Journal on optimization, 6(2):418–445, 1996.
  * 9. A. Degasperi, D. Fey, and B. N. Kholodenko.  Performance of objective functions and optimisation procedures for parameter estimation in system biology models.  NPJ systems biology and applications, 3(1):20, 2017.
  * 10. J. A. Egea, E. Balsa-Canto, M.-S. G. Garcia, and J. R. Banga.  Dynamic optimization of nonlinear processes with an enhanced scatter search method.  Industrial & Engineering Chemistry Research, 48(9):4388–4401, 2009\.
  * 11. J. A. Egea, R. Marti, and J. R. Banga.  An evolutionary method for complex-process optimization.  Computers & Operations Research, 37(2):315–324, 2010.
  * 12. A. Fiedler, S. Raeth, F. J. Theis, A. Hausser, and J. Hasenauer.  Tailored parameter optimization methods for ordinary differential equation models with steady-state constraints.  BMC systems biology, 10(1):80, 2016.
  * 13. F. Fröhlich, B. Kaltenbacher, F. J. Theis, and J. Hasenauer.  Scalable parameter estimation for genome-scale biochemical reaction networks.  PLoS computational biology, 13(1):e1005331, 2017.
  * 14. F. Fröhlich, C. Loos, and J. Hasenauer.  Scalable inference of ordinary differential equation models of biochemical processes.  In Gene Regulatory Networks, pages 385–422. Springer, New York, 2019.
  * 15. F. Fröhlich, F. J. Theis, J. O. Rädler, and J. Hasenauer.  Parameter estimation for dynamical systems with discrete events and logical operations.  Bioinformatics, 33(7):1049–1056, 2016.
  * 16. A. Gábor and J. R. Banga.  Robust and efficient parameter estimation in dynamic models of biological systems.  BMC systems biology, 9(1):74, 2015.
  * 17. H. Hass, C. Loos, E. Raimúndez-Álvarez, J. Timmer, J. Hasenauer, and C. Kreutz.  Benchmark problems for dynamic modeling of intracellular processes.  Bioinformatics, jan 2019.
  * 18. S. Hengl, C. Kreutz, J. Timmer, and T. Maiwald.  Data-based identifiability analysis of non-linear dynamical models.  Bioinformatics, 23(19):2612–2618, 2007.
  * 19. J. P. Ioannidis.  Meta-research: Why research on research matters.  PLoS biology, 16(3):e2005468, 2018.
  * 20. C. Kreutz.  New concepts for evaluating the performance of computational methods.  IFAC-PapersOnLine, 49(26):63–70, 2016.
  * 21. C. Kreutz.  An easy and efficient approach for testing identifiability.  Bioinformatics, 34(11):1913–1921, 2018.
  * 22. C. Li, M. Donizelli, N. Rodriguez, H. Dharuri, L. Endler, V. Chelliah, L. Li, E. He, A. Henry, M. I. Stefan, et al.  Biomodels database: An enhanced, curated and annotated resource for published quantitative kinetic models.  BMC systems biology, 4(1):92, 2010.
  * 23. S. Mangul, L. S. Martin, B. L. Hill, A. K.-M. Lam, M. G. Distler, A. Zelikovsky, E. Eskin, and J. Flint.  Systematic benchmarking of omics computational tools.  Nature communications, 10, 2019.
  * 24. P. Meyer, T. Cokelaer, D. Chandran, K. H. Kim, P.-R. Loh, G. Tucker, M. Lipson, B. Berger, C. Kreutz, A. Raue, et al.  Network topology and parameter estimation: from experimental design methods to gene regulatory network kinetics using a community based approach.  BMC systems biology, 8(1):13, 2014.
  * 25. C. G. Moles, P. Mendes, and J. R. Banga.  Parameter estimation in biochemical pathways: a comparison of global optimization methods.  Genome research, 13(11):2467–2474, 2003.
  * 26. M. Peifer and J. Timmer.  Parameter estimation in ordinary differential equations for biochemical processes using the method of multiple shooting.  IET Systems Biology, 1(2):78–88, 2007.
  * 27. B. Peters, S. E. Brenner, E. Wang, D. Slonim, and M. G. Kann.  Putting benchmarks in their rightful place: The heart of computational biology, 2018.
  * 28. A. Raue, C. Kreutz, T. Maiwald, J. Bachmann, M. Schilling, U. Klingmüller, and J. Timmer.  Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood.  Bioinformatics, 25(15):1923–1929, 2009.
  * 29. A. Raue, M. Schilling, J. Bachmann, A. Matteson, M. Schelke, D. Kaschek, S. Hug, C. Kreutz, B. D. Harms, F. J. Theis, et al.  Lessons learned from quantitative dynamical modeling in systems biology.  PloS one, 8(9):e74335, 2013.
  * 30. A. Raue, B. Steiert, M. Schelker, C. Kreutz, T. Maiwald, H. Hass, J. Vanlier, C. Toensing, L. Adlung, R. Engesser, W. Mader, T. Heinemann, J. Hasenauer, M. Schilling, T. Hoefer, E. Klipp, F. Theis, U. Klingmueller, B. Schoeberl, and J. Timmer.  Data2Dynamics: a modeling environment tailored to parameter estimation in dynamical systems.  Bioinformatics, 31:3558–3560, Nov. 2015.
  * 31. T. Schweiger.  The performance of constrained optimization for ordinary differential equation models.  Bachelor thesis, University of Freiburg, 2017.
  * 32. B. Steiert, A. Raue, J. Timmer, and C. Kreutz.  Experimental design for parameter estimation of gene regulatory networks.  PloS one, 7(7):e40052, 2012.
  * 33. J. Sun, J. M. Garibaldi, and C. Hodgman.  Parameter estimation using metaheuristics in systems biology: a comprehensive review.  IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB), 9(1):185–202, 2012.
  * 34. A. F. Villaverde, F. Fröhlich, D. Weindl, J. Hasenauer, and J. R. Banga.  Benchmarking optimization methods for parameter estimation in large kinetic models.  Bioinformatics, 35(5):830–838, 2018.
  * 35. A. F. Villaverde, D. Henriques, K. Smallbone, S. Bongard, J. Schmid, D. Cicin-Sain, A. Crombach, J. Saez-Rodriguez, K. Mauch, E. Balsa-Canto, P. Mendes, J. Jaeger, and J. R. Banga.  BioPreDyn-bench: a suite of benchmark problems for dynamic modelling in systems biology.  BMC Syst. Biol., 9:8, 2015.
  * 36. L. M. Weber, W. Saelens, R. Cannoodt, C. Soneson, A. Hapfelmeier, P. P. Gardner, A.-L. Boulesteix, Y. Saeys, and M. D. Robinson.  Essential guidelines for computational method benchmarking.  Genome Biology, 20(1):125, 2019.
  * 37. F. Wieland.  Implementation and assessment of optimization approaches for parameter estimation in systems biology.  Masters thesis, University of Freiburg, 2016.
  * 38. B. Winkelmann.  Implementatoin and evaluation of optimization algorithms for parameter estimation in systems biology.  Bachelor thesis, University of Freiburg, 2016.


# Lipschitz continuity for harmonic functions and solutions of the
$\bar{\alpha}$-Poisson equation

Miodrag Mateljević, Nikola Mutavdžić, Adel Khalfallah  Faculty of mathematics,
University of Belgrade, Studentski Trg 16, 11000 Belgrade, Serbia
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We study Lipschitz continuity for solutions of the $\bar{\alpha}$-Poisson
equation in planar cases. We also review some recently obtained results. As
corolary we can restate results for harmonic and gradient harmonic functions.

11footnotetext:  Mathematics Subject Classification 2010 Primary 42B15,
Secondary 42B30. 22footnotetext:  Key words and phrases: Poisson kernel,
harmonic functions.

##  1\. Introduction and preliminaries

The weighted Laplacian operator $L_{\rho}$ is defined by $L_{\rho}=D_{z}(\rho
D_{\overline{z}})$ and $L^{*}_{\rho}=D_{\overline{z}}(\rho D_{z})$. If the
weight function $\rho$ is chosen to be
$\rho_{\alpha}(z)=(1-|z|^{2})^{-\alpha}$ ($\alpha>-1$) in the unit disk
$\mathbb{U}$, we call $L_{\rho}$ the standard weighted Laplacian operator and
write it by $L_{\alpha}$ for simplicity, and $\overline{L_{\alpha}}$ will be
notation for $L^{*}_{\rho}$.

Our main result is:

######  Theorem 8.

Assume that $g\in C(\mathbb{D})$ is such that $(1-|z|^{2})^{\alpha}g$ is
bounded and $u\in V_{\mathbb{D}\rightarrow\Omega}[g]$ with the representation
$u(w)=v(w)+G_{\alpha}[g](w)$. If $\alpha>0$, and $\underline{v}$ is Lipshitz,
then $u$ is also Lipshitz continuous on ${\mathbb{D}}$.

We can restate this result by means of certain solutions to $\alpha-$Poisson’s
equation. First we consider some basic properties of $\alpha-$harmonic
mappings. In particular, we improve Chen and Kalaj result [6]. Behm [2] found
Green function and solved the Dirichlet boundary value problem of the
$\alpha$-Poisson equation. Our method is based on Theorem 7 which gives
estimate of the Green potential $G_{\alpha}$ of $g$ and the local
$C^{2}$-coordinate method flattering the boundary [23].

At the begining of this paper we will introduce basic notation together with
definition of so called $\alpha-$Laplacian and $\alpha-$harmonic functions.
Also, definition and properties of $\alpha-$Poisson’s kernel and
$\alpha-$Poisson’s integral are stated, as a very important technical assets
used in our research. More information about this notion can be found in
Olofsson’s and Wittsten’s paper [29]. In the proceeding we recall definition
of Green function for $\alpha-$Laplacian which is thoroughly invesigated in
Behm [2]. Formulation and solution of Dirichlet boundary value problem for
$\alpha-$Poisson’s equation, proven in Chen and Kalaj paper [3], is shown in
Theorem 1. In paper [7] Chen used this result to prove the boundary
characterizations of a Lipschitz continuous $\overline{\alpha}$-harmonic
mappings, and proved Theorem 2.

As an introductory result of this paper, we loose assumption on boundary value
of $\overline{\alpha}$-harmonic mapping $v$, which is written in part $(d)$ of
Theorem 2 and attain Theorem 5. Proof of this theorem is based on Hardy space
technique which can be found in the first author’s monography [27], and
Theorem 6 proven in first author’s and A. Khalfallah’s paper [13]. The second
improvement of Theorem 2 consider the condition on
$g=-\overline{L_{\alpha}}u$. This result is proven in Theorem 7, and uses
various estimates which we establish in Sunsection 2.3.

###  1.1. $\alpha-$harmonic mappings

Linear and semilinear equations can be treated together. We take
$a(x,y)u_{x}+b(x,y)u_{y}=c(x,y,u)$, where a, b, c are $C^{1}$ functions of
their arguments. The operator $a(x,y)u_{x}+b(x,y)u_{y}$ on the left hand side
of this equation represents differentiation in a direction ($a,b$) at the
point $(x,y)$ in $(x,y)$-plane. Let us consider a curve, whose tangent at each
point ($x,y$) has the direction $(a,b)$. Coordinates $(x(s),y(s))$ of a point
on this curve satisfy (1) $dx/ds=a(x,y),dy/ds=b(x,y)$ or (2)
$dy/dx=b(x,y)/a(x,y)$.

Two complex derivatives $\frac{\partial}{\partial z}=D_{z}$ and
$\frac{\partial}{\partial\overline{z}}=D_{\overline{z}}$ of $u$ are written by

$\frac{\partial}{\partial
z}u=D_{z}u=\frac{1}{2}(u_{x}-iu_{y})\quad\mbox{and}\quad\frac{\partial}{\partial\overline{z}}u=D_{\overline{z}}u=\frac{1}{2}(u_{x}+iu_{y})$

respectively, where $z=x+iy$.

The weighted Laplacian operator $L_{\rho}$ is defined by

$L_{\rho}=D_{z}(\rho D_{\overline{z}})\quad\mbox{and}\quad
L^{*}_{\rho}=D_{\overline{z}}(\rho D_{z}).$

If the weight function $\rho$ is chosen to be
$\rho_{\alpha}(z)=(1-|z|^{2})^{-\alpha}$ in the unit disk $\mathbb{U}$, we
call $L_{\rho}$ the standard weighted Laplacian operator and write it by
$L_{\alpha}$ for simplicity, here $\alpha$ is a real number with $\alpha>-1$.
It is clear that

$\displaystyle L^{*}_{\rho}u=\overline{L_{\alpha}}u=D_{\overline{z}}\rho
D_{z}u+\rho
D_{\overline{z}z}u=\alpha(1-|z|^{2})^{-\alpha-1}zu_{z}+(1-|z|^{2})^{-\alpha}u_{\overline{z}z}$
$\displaystyle
L_{\alpha}u=\alpha(1-|z|^{2})^{-\alpha-1}\overline{z}u_{\overline{z}}+(1-|z|^{2})^{-\alpha}u_{\overline{z}z}$

First we can see that $L^{*}_{\alpha}u=0$ iff

$(1-|z|^{2})u_{\overline{z}z}+\alpha zu_{z}=0$ (1)

It can be easily checked that $L_{\alpha}u=0$ iff
$L^{*}_{\alpha}\overline{u}=0.$

Set $p=u_{z}$ and $q=u_{\overline{z}}$. Since
$\overline{u}_{\overline{z}}=\overline{u_{z}}$ and
$\overline{u}_{z}=\overline{u_{\overline{z}}}$, we find

$zu_{z}$ and $\overline{z}u_{\overline{z}}$ are conjugate, and also

$u_{\overline{z}z}=q_{z}$ and
$\overline{u}_{\overline{z}z}=\overline{q}_{\overline{z}}=\overline{q_{z}}$
and therefore $u_{\overline{z}z}$ and $\overline{u}_{\overline{z}z}$ are
conjugate.

If we set $d(z)=1-|z|^{2}$, then $\rho_{\alpha}=d^{-\alpha}$ and by easy
computation we find

$\rho_{z}=\alpha d^{-\alpha}\overline{z},\quad\rho_{\overline{z}}=\alpha
d^{-\alpha-1}z,\quad\rho_{x}=2\alpha
d^{-\alpha-1}x\quad\mbox{and}\quad\rho_{y}=2\alpha d^{-\alpha-1}y.$

Since, $2\rho D_{z}u=\rho(u_{x}-iu_{y})$, we find

$4L_{\rho}=D_{x}[\rho(u_{x}-iu_{y})]+iD_{y}[\rho(u_{x}-iu_{y})]=D_{x}(\rho
u_{x})+D_{y}(\rho u_{y})+i(\rho_{y}u_{x}-\rho_{x}u_{y}).$

Hence

$4L_{\rho}=\rho\Delta
u+\rho_{x}u_{x}+\rho_{y}u_{y}+i(\rho_{y}u_{x}-\rho_{x}u_{y}).$

If $u$ is real-valued function then $L_{\rho}u=0$ iff $\rho\Delta
u+\rho_{x}u_{x}+\rho_{y}u_{y}=0$ and $yu_{x}-xu_{y}=0$, that is

$\Delta u+2\alpha\rho_{1}(xu_{x}+yu_{y})\mbox{ and }yu_{x}-xu_{y}=0.$

By [38], solutions of $yu_{x}-xu_{y}=0$ is $u=f(x^{2}+y^{2})$. Since $\rho
u_{z}=\rho g(r)\overline{z}$, we find $\rho g(r)r^{2}=zF(z)=c$ and hence $F=0$
and $u_{z}=0$. Thus $u=c$.

If a function $u\in C^{2}(\mathbb{U})$ satisfies the $\alpha$-harmonic
equation

$L_{\alpha}(u)=0,$

then we call it an $\alpha$-harmonic mapping. In the case $\alpha=0$,
$\alpha$-harmonic mappings are just Euclidean harmonic mappings. For
$L_{\alpha}$ notation $\Delta_{\alpha}$ is also used in the literature.

###  1.2. $\alpha$-Poisson’s integral

In [29], Olofsson and Wittsten showed that if an $\alpha$-harmonic function
$f$ satisfies

$\lim_{r\to
1^{-}}f_{r}=f^{\ast}\in\mathcal{D}^{\prime}(\mathbb{T})\;\;(\alpha>-1),$

then it has the form of a _Poisson type integral_

$\displaystyle
f(z)=P_{\alpha}[f^{\ast}](z)=\frac{1}{2\pi}\int^{2\pi}_{0}P_{\alpha}(ze^{-i\theta})f^{\ast}(e^{i\theta})d\theta$

in $\mathbb{D}$, where

$P_{\alpha}(z)=\frac{(1-|z|^{2})^{\alpha+1}}{(1-z)(1-\overline{z})^{\alpha+1}}$

is the $\alpha$-harmonic (complex valued) Poisson kernel in ${\mathbb{D}}$. In
the case $\alpha=0$ we have classiacal Poisson’s kernel for harmonic function
and we write it as $P$ instead of $P_{0}$. Also, we write $P[f^{*}]$ instead
of $P_{0}[f^{*}]$ for Piosson’s integral of the function $f^{*}$.

##  2\. Lipschitz continuity for solutions of the $\bar{\alpha}$-Poisson
equation

###  2.1. An introductory result

As a starting point of our investigation we used Theorem 2 which can be found
in Chen’s paper [7]. This theroem gives some rather strong assumption on
$g=-\overline{L_{\alpha}}u$ ($g\in C({\mathbb{D}}))$ as well as for boundary
values of $u$ (condition $(d)$ of Theorem 2), which are proven to be
sufficient for Lipshitz continuity of $u$.

In order to formulate basic result we need to introduce some preliminary
notes. Let $V_{D\rightarrow\Omega}[g]$ denote the family of solutions of
$v:D\rightarrow\Omega$, $v\in C^{2}(\mathbb{D})$ of the
$\overline{\alpha}$-Poisson equation

$\left\\{\begin{array}[]{ll}v(z)=f(z),&\hbox{if}\,\,z\in\mathbb{T},\\\
-(\overline{L_{\alpha}})v(z)=g,&\hbox{if}\,\,z\in\mathbb{D},\end{array}\right.$
(2)

where $g\in C(\overline{\mathbb{D}})$, $f\in L^{1}(\mathbb{D})$ is the limit
of $v(re^{i\theta})$ as r tends to $1^{-}$, and $v$ is a sense-preserving
diffeomorphism.

For the case wherein the boundary function $f$ vanishes, Behm [2] solved the
above Dirichlet boundary value problem of the $\overline{\alpha}-$Poisson
equation. In paper [3], Chen and Kalaj combined the representation theorem
given by Olofsson and Wittsten [29] with the one given by Behm. They obtained
the following theorem.

######  Theorem 1 ([3]).

Let $g\in C(\mathbb{D})$ be such that $(1-|z|^{2})^{\alpha+1}g$ belongs to
$L^{1}(\mathbb{D})$, wherein $\alpha>-1$. If $u\in C^{2}(\mathbb{D})$ is a
solution of the equation $-\overline{L_{\alpha}}u=g$ satisfying the condition
$u(re^{it})$ converges to a function $f\in L^{1}(\mathbb{T})$ as $r\to 1^{-}$,
then for every $w\in\mathbb{D}$ we have

$u(w)=v(w)+G_{\alpha}[g](w),\vspace{-2mm}$

where

$v(w){=}\frac{1}{2\pi}\int_{\mathbb{T}}\frac{(1{-}|w|^{2})^{\alpha+1}}{(1{-}z\overline{w})(1{-}\overline{z}w)^{\alpha+1}}f(z)\,\mathrm{d}\theta,\quad
G_{\alpha}[g](w){=}\iint_{\mathbb{D}}G_{\alpha}(z,w)g(z)\,\mathrm{d}x\,\mathrm{d}y,$

and the Green function $G_{\alpha}(z,w)$ of the adjoint Laplacian $L_{\alpha}$
is given by

$\displaystyle
G_{\alpha}(z,w)=\frac{(1-\overline{z}w)^{\alpha}h(q(z,w))}{2\pi},\mbox{ with
}z\neq w,$ $\displaystyle
h(r)=\frac{1}{2}\int_{0}^{1-r^{2}}\frac{t^{\alpha}}{1-t}\,\mathrm{d}t,q(z,w)=\left|\frac{z-w}{1-\overline{w}z}\right|.$

In [7] Chen provided the boundary characterizations of a Lipschitz continuous
$\overline{\alpha}$-harmonic mapping as follows

Define $\underline{f}(t)=f(e^{it})$ and

$S[f](w)=\frac{1}{\pi}\int_{0}^{2\pi}\frac{(1-|w|^{2})^{\alpha}}{(1-z\overline{w})^{\alpha}}\frac{Im(w\overline{z})}{|z-w|^{2}}\underline{f}(t)dt$
(3)

where $z=e^{it}$.

######  Theorem 2 ([7]).

Assume that $g\in C(\overline{\mathbb{D}})$ and $u\in
V_{\mathbb{D}\rightarrow\Omega}[g]$ with the representation
$u(w)=v(w)+G_{\alpha}[g](w)$. If $\alpha\geqslant 0$, then the following
conditions are equivalent:

  * $(a)$

$u$ is a (K, K’)-quasiconformal mapping and $|\frac{\partial}{\partial
r}v|\leq L$ on $\mathbb{D}$, and $L$ is a constant.

  * $(b)$

$u$ is Lipschitz continuous with the Euclidean metric.

  * $(c)$

$v$ is Lipschitz continuous with the Euclidean metric.

  * $(d)$

$f$ is absolutely continuous on $\mathbb{T}$, $f^{\prime}\in L^{\infty}$ and
$S[\underline{f}^{\prime}]$ is bounded on $\mathbb{D}$.

In order to prove the main result of this paper we will need to prove two
refinement of the result above.

###  2.2. Refinement of part $(d)$ in Theorem 2

As we can see in [13], for $p\in(0,\infty]$, the generalized Hardy space
$\mathcal{H}^{p}_{\mathcal{G}}(\mathbb{D})$ consists of all measurable
functions from $\mathbb{D}$ to $\mathbb{C}$ such that $M_{p}(r,f)$ exists for
all $r\in(0,1)$, and $\|f\|_{p}<\infty$, where

$M_{p}(r,f)=\left(\frac{1}{2\pi}\int_{0}^{2\pi}|f(re^{i\theta})|^{p}\,d\theta\right)^{\frac{1}{p}}$

and

$\|f\|_{p}=\begin{cases}\displaystyle\sup\\{M_{p}(r,f):\;0<r<1\\}&\mbox{if
}p\in(0,\infty),\\\ \displaystyle\sup\\{|f(z)|:\;z\in\mathbb{D}\\}&\mbox{if
}p=\infty.\end{cases}$

The classical Hardy space $\mathcal{H}^{p}(\mathbb{D})$ (resp.
$h^{p}(\mathbb{D})$) is the set of all elements of
$\mathcal{H}^{p}_{\mathcal{G}}(\mathbb{D})$ which are analytic (resp.
harmonic) on ${\mathbb{D}}$.

######  Definition 3 ([27]).

Let $\psi\in L^{1}([0,2\pi])$. Then

$H(\psi)(\varphi)=-\frac{1}{2\pi}\int\limits_{0_{+}}^{\pi}\frac{\psi(\varphi+t)-\psi(\varphi-t)}{\tan
t/2}\,\mathrm{d}t,$

denotes the Hilbert transformations of $\psi.$

Recall that $\underline{h}(t)=h(e^{it})$. The following property of the
Hilbert transform is also sometimes taken as the definition:

######  Theorem 4 ([27]).

If $u=P[\psi]$ and $v$ is the harmonic conjugate of $u$, then
$\underline{v}(t)=H(\psi)(t)$ a.e.

Note that, if $\psi$ is $2\pi$-periodic, absolutely continuous on $[0,2\pi]$
(and therefore $\psi^{\prime}\in L^{1}[0,2\pi]$), then

$\frac{\partial h}{\partial\theta}=P[\psi^{\prime}].$

Hence, since $r\frac{\partial}{\partial r}h$ is the harmonic conjugate of
$\frac{\partial}{\partial\theta}h$, we find

$\displaystyle r\frac{\partial h}{\partial r}$
$\displaystyle=P[H(\psi^{\prime})],\quad$ (4)
$\displaystyle\left(\frac{\partial h}{\partial r}\right)^{*}(e^{i\theta})$
$\displaystyle=H(\psi^{\prime})(\theta)\quad{\rm a.e.}\,$ (5)

Using the above outline we can derive:

######  Theorem 5.

Let $h$ be harmonic on $\mathbb{U}$ and $h\in h^{1}(\mathbb{U})$. Then $h$ is
Lip iff $\underline{h}^{\prime}\in L^{\infty}$ and
$H(\underline{h}^{\prime})\in L^{\infty}$.

Let $h$ be gradient harmonic or $\alpha$-harmonic, $\alpha>0$, on $\mathbb{U}$
and $h\in h^{1}(\mathbb{U})$. Then $h$ is Lip iff $\underline{h}^{\prime}\in
L^{\infty}$.

The second part of the last theorem is a direct consequence of the first
author’s result, together with A. Khalfallah, which is stated below.

######  Theorem 6 ([13]).

Let $\alpha\in(-1,\infty)$ with $\alpha\not=0$ and let
$f=\mathcal{P}_{\alpha}[F]$ and $F$ is absolute continuous such that
$\dot{F}\in L^{p}$ with $1\leqslant p\leqslant\infty$.

  1. (1)

If $\alpha>0$, then $\frac{\partial}{\partial{\overline{z}}}f$ and
$\frac{\partial}{\partial z}f$ are in
$\mathcal{H}_{\mathcal{G}}^{p}({\mathbb{D}})\subset L^{p}({\mathbb{D}})$.

  2. (2)

If $\alpha\in(-1,0)$, then $\frac{\partial}{\partial{\overline{z}}}f$ and
$\frac{\partial}{\partial z}f$ are in $L^{p}({\mathbb{D}})$ for
$p<-\frac{1}{\alpha}$.

  3. (3)

For $\alpha\in(-1,0)$ and $p\geqslant-\frac{1}{\alpha}$ there exists $f$ an
$\alpha$-harmonic function such that
$\frac{\partial}{\partial{\overline{z}}}f$ and $\frac{\partial}{\partial z}f$
$\not\in L^{p}({\mathbb{D}})$. Moreover,
$\frac{\partial}{\partial{\overline{z}}}f$ and $\frac{\partial}{\partial z}f$
$\not\in\mathcal{H}_{\mathcal{G}}^{1}({\mathbb{D}})$.

###  2.3. Refinement of the condition on $g=-\overline{L_{\alpha}}u$ in
Theorem 2

In this subsection we will prove that instead of $g\in
C(\overline{\mathbb{D}})$ we can use assume that $g\in C(\mathbb{D})$ can be
such that $|g(z)|\leqslant M(1-|z|^{2})^{-\alpha},z\in\mathbb{D}$, in oreder
to prove Lipshitz continuity of $\overline{\alpha}$-Green integral
$G_{\alpha}[g]$ of function $g$. This fact will play an important part in the
proof of our main result.

The following two estimates can be obtained by direct investigation of the
Green function $G_{\alpha}$, and can be found in [7].

$\displaystyle 2\pi\left|\frac{\partial}{\partial w}G_{\alpha}(z,w)\right|$
$\displaystyle\leqslant\alpha
C_{\alpha}|1-\overline{z}w|^{\alpha-1}\left(1{-}\left|\frac{z{-}w}{1{-}\overline{w}z}\right|^{2}\right)^{\alpha{+}1}\left(1{-}\log\left|\frac{z{-}w}{1{-}\overline{w}z}\right|^{2}\right)$
$\displaystyle+\frac{(1-|z|^{2})^{\alpha+1}(1-|w|^{2})^{\alpha}}{2|1-z\overline{w}|^{\alpha+1}|z-w|},\vspace{-2mm}$
(6) $\displaystyle
2\pi\left|\frac{\partial}{\partial\overline{w}}G_{\alpha}(z,w)\right|\leqslant\frac{(1-|z|^{2})^{\alpha+1}(1-|w|^{2})^{\alpha}}{2|1-\overline{z}w|^{\alpha+1}|z-w|}.$
(7)

In order to start with our work, we will prove the following two tehnical
lemmas.

######  Lemma 1.

If $\beta>1$ then

$\int_{0}^{2\pi}\frac{\,\mathrm{d}t}{|1-r\rho
e^{it}|^{\beta}}\preceq\frac{1}{|1-r\rho|^{\beta-1}}$

for $0<r,\rho<1$.

###### Proof.

$\displaystyle\int_{0}^{2\pi}\frac{\,\mathrm{d}t}{|1-r\rho e^{it}|^{\beta}}$
$\displaystyle=2\int_{0}^{\pi}\frac{\,\mathrm{d}t}{((1-r\rho)^{2}+4r\rho\sin^{2}\frac{t}{2})^{\beta/2}}$
$\displaystyle\leqslant\int_{0}^{\pi}\frac{\,\mathrm{d}t}{((1-r\rho)^{2}+c_{1}t^{2})^{\beta/2}}\leqslant\left|t=(1-r\rho)u\right|$
$\displaystyle\leqslant\int_{0}^{\pi/(1-r\rho)}\frac{(1-r\rho)\,\mathrm{d}u}{(1-r\rho)^{\beta}(1+c_{1}u^{2})^{\beta/2}}$
$\displaystyle\leqslant\frac{1}{(1-r\rho)^{\beta-1}}\int_{0}^{\infty}\frac{\,\mathrm{d}u}{(1+c_{1}u^{2})^{\beta/2}}.$

since the last integral converges, we have desired result. ∎

######  Lemma 2.

There exists $c_{2}>0$ such that

$M_{1}(r)=\iint_{\mathbb{D}}\frac{\,\mathrm{d}x\,\mathrm{d}y}{|z-r|}\leqslant
c_{2}$

for every $0<r<1$.

###### Proof.

Let us use the supsitution $z-r=\rho e^{it}$, where $0\leqslant
t<2\pi,0<\rho<\rho(t)=|r-e^{it}|\leqslant r+1$. Then

$\displaystyle\iint_{\mathbb{D}}\frac{\,\mathrm{d}x\,\mathrm{d}y}{|z-r|}=\int_{0}^{2\pi}\,\mathrm{d}t\int_{0}^{\rho(t)}\frac{\rho\,\mathrm{d}\rho}{\rho}\leqslant\int_{0}^{2\pi}(r+1)\,\mathrm{d}t\leqslant
4\pi.$ (8)

∎

Let $|w|=r$,

$\displaystyle I_{1}(w)$
$\displaystyle=\iint_{\mathbb{D}}\frac{(1-|z|^{2})^{\alpha+1}(1-|w|^{2})^{\alpha}}{2|1-\overline{z}w|^{\alpha+1}|z-w|}\,\mathrm{d}x\,\mathrm{d}y,$
$\displaystyle I_{2}(w)$
$\displaystyle=\iint_{\mathbb{D}}|1{-}\overline{z}w|^{\alpha-1}\left(1{-}\left|\frac{z{-}w}{1{-}\overline{w}z}\right|^{2}\right)^{\alpha{+}1}\left(1{-}\log\left|\frac{z{-}w}{1{-}\overline{w}z}\right|^{2}\right)\,\mathrm{d}x\,\mathrm{d}y.$

Also, inequalities

$|1-\overline{w}\zeta|\geqslant
1-|w|\quad\mbox{and}\quad|1-\overline{w}\zeta|\geqslant 1-|\zeta|$ (9)

can easily be veryfied.

The following two lemmas are crucial for the main result of this section:

######  Lemma 3.

There exists $c_{3}>0$ such that

$I_{1}(w)\leqslant c_{3}(1-|w|^{2})^{\alpha}$ (10)

for every $|w|<1$.

###### Proof.

Using (9) we get

$I_{1}(w)\preceq(1-|w|^{2})^{\alpha}\iint_{\mathbb{D}}\frac{\,\mathrm{d}x\,\mathrm{d}y}{|z-w|}.$

Since we can use coordinate change $s=\frac{w}{|w|}z$, we can use Lemma 2 to
get our result. ∎

Let $\zeta=\varphi_{w}(z)=\frac{w-z}{1-\overline{w}z}$. Then, for each
$w\in\mathbb{D},\varphi_{w}$ is a conformal mapping of $\mathbb{D}$ satisfying
the following identities:

$\displaystyle z$
$\displaystyle=\frac{w-\zeta}{1-\overline{w}\zeta},1-|z|^{2}=\frac{(1-|z|^{2})(1-|w|^{2})}{|1-\overline{w}\zeta|^{2}},$
(11) $\displaystyle 1$
$\displaystyle-\overline{w}z=\frac{1-|w|^{2}}{1-\overline{w}\zeta},\,\mathrm{d}z=-\frac{1-|w|^{2}}{(1-\overline{w}\zeta)^{2}}\,\mathrm{d}\zeta$

######  Lemma 4.

There exists $c_{4}>0$ such that

$I_{2}(w)\leqslant c_{4}(1-|w|^{2})^{\alpha}$ (12)

for every $|w|<1$.

###### Proof.

By using supstitution $s=\frac{w}{|w|}\zeta$, and $s=\rho e^{it}$ we get

$\displaystyle I_{2}(w)$
$\displaystyle=I_{2}(r)=\iint_{\mathbb{D}}\frac{(1-|w|^{2})^{\alpha+1}(1-|\zeta|^{2})^{\alpha+1}}{|1-\overline{w}\zeta|^{\alpha+3}}(1-\log|\zeta|^{2})\,\mathrm{d}\xi\,\mathrm{d}\eta$
$\displaystyle=(1{-}|w|^{2})^{\alpha}\int_{0}^{1}(1{-}\rho^{2})^{\alpha+1}(1{-}r^{2})(1{-}\log\rho^{2})\int_{0}^{2\pi}\frac{\,\mathrm{d}t}{|1{-}r\rho
e^{it}|^{\alpha+3}}\rho\,\mathrm{d}\rho.$

Using Lemma 1 we get

$\displaystyle I_{2}(r)$
$\displaystyle\preceq(1-|w|^{2})^{\alpha}\int_{0}^{1}\frac{(1-\rho^{2})^{\alpha+1}(1-r^{2})(1-\log\rho^{2})}{|1-r\rho|^{\alpha+2}}\rho\,\mathrm{d}\rho.$

Since $1-r\rho\geqslant 1-r$ and $1-r\rho\geqslant 1-\rho$ we have that

$\displaystyle
I_{2}(r)\preceq(1-|w|^{2})^{\alpha}\int_{0}^{1}\rho(1-\log\rho^{2})\,\mathrm{d}\rho\leqslant
c_{4}(1-|w|^{2})^{\alpha}$

for some $c_{4}>0$ which does not depent of $0\leqslant r<1$. ∎

We are now ready to formulate the main result of this section, which is
generalisation of Lemma 3.4 in Chen’s paper [7]. The proof of this result
follows directly from Lemma 3 and Lemma 4.

######  Theorem 7.

Let $g\in C(\mathbb{D})$ be such that $|g(z)|\leqslant
M(1-|z|^{2})^{-\alpha},z\in\mathbb{D}$ for some $M>0$ and let $\alpha>0$ be
arbitrary. Assume that $G_{\alpha}[g]$ is the Green potential of $g$ given by

$G_{\alpha}[g](w)=\iint_{\mathbb{D}}G_{\alpha}(z,w)g(z)\,\mathrm{d}x\,\mathrm{d}y.$

Then $\frac{\partial}{\partial w}G_{\alpha}[g]$ and
$\frac{\partial}{\partial\overline{w}}G_{\alpha}[g]$ are both bounded in the
unit disc $\mathbb{D}$.

As a ditect consequence of Theorem 7 and Theorem 5 we have the main result of
this paper.

######  Theorem 8.

Assume that $g\in C(\mathbb{D})$ is such that $(1-|z|^{2})^{\alpha}g$ is
bounded and $u\in V_{\mathbb{D}\rightarrow\Omega}[g]$ with the representation
$u(w)=v(w)+G_{\alpha}[g](w)$. If $\alpha>0$, and $\underline{v}$ is Lipshitz,
then $u$ is also Lipshitz continuous on ${\mathbb{D}}$.

## References

  * [1] M. Arsenović, M. Mateljević, Modulus of continuity of normal derivative of a harmonic functions at a boundary point preprint 2022, accepted in Filomat, March 2023.
  * [2] Behm G. Solving Poisson equation for the standard weighted Laplacian in the unit disc. ArXiv:1306.2199v2, 2013.
  * [3] X. D. Chen, D. Kalaj, Representation theorem of standard weighted harmonic mappings in the unit disk. Https:// www.researchgate.net/publication/330243622, 2018.
  * [4] A. Gjokaj, D. Kalaj, _QCH mappings between unit ball and domain with $C^{1,\alpha}$ boundary_, 2021, arXiv:2005.05667 [math.AP].
  * [5] A. Dall’Acqua, G. Sweers, _Estimates for Green function and Poisson kernels of higher order Dirichlet boundary value problems,_ Journal of Differential Equations, 205, no. 2, 2004, p. 466–487.
  * [6] X. Chen, D. Kalaj, _A representation theorem for standard weighted harmonic mappings with an integer exponent and its applications_ , J. Math. Anal. Appl. 444, 2 (2016), 1233–1241.
  * [7] X. Chen, _Lipschitz continuity for solutions of the $\bar{\alpha}$-Poisson equation_, Sci. China, Math. 62, 10 (2019), 1935–1946.
  * [8] J. Chen, M. Huang, A. Rasila, X. Wang, _On Lipschitz continuity of solutions of hyperbolic Poisson’s equation,_ Calc. Var. Partial Differ. Equ., 57, no. 1, 2018, p. 1–32.
  * [9] S. L. Chen and S. Ponnusamy, Schwarz’s lemmas for mappings satisfying Poisson’s equation, Indag. Math., New Ser., 30 (2019), 1087-1098.
  * [10] S. L. Chen, S. Ponnusamy, X. T. Wang,Remarks on ’Norm estimates of the partial derivatives for harmonic mappings and harmonic quasiregular mappings, J. Geom. Anal. 31, 11051-11060 (2021).
  * [11] S. L. Chen, The Heinz type inequality, Bloch type theorem and Lipschitz characteristic of polyharmonic mappings, to appear, March 2023.
  * [12] D. Gilbarg and N. Trudinger, Elliptic partial Differential Equation of Second Order, Springer-Verlag, Berlin, Second Edition, 1983.
  * [13] A. Khalfallah, M. Mateljević, Norm estimates of the partial derivatives and Schwarz lemma for $\alpha$-harmonic functions, Accepted in Complex Variables and Elliptic Equations.
  * [14] B. Hollenbeck and I. E. Verbitsky, Best Constants for the Riesz Projection, J. Fun. Anal., 175 (2), 370-392. (2000)
  * [15] D. Kalaj, On Riesz type inequalities for harmonic mappings on the unit disk, Trans. Am. Math. Soc. 372, No. 6, 4031-4051 (2019).
  * [16] P. Li, L. Tam, Uniqueness and Regularity of Proper Harmonic Maps II, Indiana University Mathematics Journal, 42 (2), 591–635 (1993).
  * [17] M. Mateljević,  The Lower Bound for the Modulus of the Derivatives and Jacobian of Harmonic Injective Mappings. \- Filomat 29:2, 2015, 221-244.
  * [18] M. Mateljević, V. Božin, M. Knežević: Quasiconformality of harmonic mappings between Jordan domains, Filomat, Vol 24, No 3, 2010, 111-124.
  * [19] M. Mateljević, N. Mutavdžić, On Lipschitz continuity and smoothness up to the boundary of solutions of hyperbolic Poisson’s equation, arXiv:2208.06197v1 [math.CV] 12 Aug 2022.
  * [20] L. Ma, _Hölder continuity of hyperbolic Poisson integral and hyperbolic Green integral_ , Monatshefte für Mathematik, 199, 2022,
  * [21] M. Mateljević, N. Mutavdžić, On Lipschitz continuity and smoothness up to the boundary of solutions of hyperbolic Poisson’s equation, submitted for publication, September 2022.
  * [22] M. Mateljević, N. Mutavdžić, The Boundary Schwarz lemma for harmonic and pluriharmonic mappings and some generalizations, accepted in Bulletin of the Malaysian Mathematical Sciences Society, June 2022.
  * [23] M. Mateljević, _Boundary behaviour of partial derivatives for solutions to certain Laplacian-gradient inequalities and spatial qc maps_ , Springer Proceedings in Mathematics & Statistics, 357, 2021, p. 393-418.
  * [24] M. Mateljević, _Boundary Behaviour Of Partial Derivatives For Solutions To Certain Laplacian-Gradient Inequalities And Spatial Qc Maps 2,_ Preprint, Comunicated at XII Symposium Mathematics and Applications, Mathematical Faculty, Belgrade, 2,3, 2022.
  * [25] M. Mateljević, R. Salimov, E. Sevost’yanov, _Hölder and Lipschitz continuity in Orlicz-Sobolev classes, distortion and harmonic mappings,_ Filomat, 32 no. 16 p. 5361–5392.
  * [26] M. Mateljević, E. Sevost’yanov, On the behavior of Orlicz-Sobolev mappings with branching on the unit sphere, Ukr. Math. Bull 19(2022), No 4, 542-584; Journal of Mathematical Sciences volume 270, pages 467-499 (2023).
  * [27] M. Mateljević, Topics in Conformal, Quasiconformal and Harmonic maps, Zavod za udzbenike, Beograd, 2012.
  * [28] C. A. Nodler, D. M. Oberlin, Moduli of continuity and a Hardy-Littlewood theorem, Lecture Notes in Math. 1351, p. 265-272, Springer-Verlag, Berlin etc., 1988.
  * [29] A. Olofsson and J. Wittsten, Poisson integrals for standard weighted Laplacians in the unit disc J. Math. Soc. Japan, 65 (2013), 447–486.
  * [30] M. Riesz, Sur les fonctions conjugées, Math. Zeit. 27, 218-244 (1927).
  * [31] W. Rudin, Real and complex analysis, McGraw-Hill Book Co, 1966.
  * [32] W. Rudin, Function Theory in the Unit Ball, Springer-Verlag, New York, 1980.
  * [33] E. M. Stein, Singular Integrals and Differentiability Properties of Functions,, Princeton University Press, 1970.
  * [34] E. M. Stein, G. Weiss. Introduction to Fourier analysis on Euclidean spaces, Princeton University Press, 1971.
  * [35] M. Stoll, Harmonic and Subharmonic Function Theory on the Hyperbolic Ball (London Mathematical Society Lecture Note Series), Cambridge: Cambridge University Press (2016).
  * [36] J. F. Zhu, Norm estimates of the partial derivatives for harmonic mappings and harmonic quasiregular mappings, J. Geom. Anal., 31, 5505-5525 (2021)
  * [37] W. Kjell-Ove, _Inequalities for the Green Function and Boundary Continuity of the Gradient of Solutions of Elliptic Differential Equations_ , Mathematica Scandinavica, 21, No. 1 , 1967, p. 17–37.
  * [38] https://en.wikipedia.org/wiki/First-order_partial_differential_equation